Draft version October 28, 2010 

Preprint typeset using LAT^X style cmulatcapj v. 11/10/09 



SATURATED TORQUE FORMULA FOR PLANETARY MIGRATION IN VISCOUS DISKS WITH THERMAL 
DIFFUSION: RECIPE FOR PROTOPLANET POPULATION SYNTHESIS 

F. S. Masset 1 ' 2 

Laboratoire AIM, CEA/DSM - CNRS - Universite Paris Diderot, Irfu/Service d'Astrophysique, Bat. 709, CEA/Saclay, 91191 

Gif-sur-Yvette, France 

AND 

J. Casoli 

Laboratoire AIM, CEA/DSM - CNRS - Universite Paris Diderot, Irfu/Service d'Astrophysique, Bat. 709, CEA/Saclay, 91191 

Gif-sur-Yvette, France 
Draft version October 28, 2010 

ABSTRACT 

We provide torque formulae for low mass planets undergoing type I migration in gaseous disks. 
These torque formulae put special emphasis on the horseshoe drag, which is prone to saturation: the 
asymptotic value reached by the horseshoe drag depends on a balance between coorbital dynamics 
(which tends to cancel out or saturate the torque) and diffusive processes (which tend to restore 
the unperturbed disk profiles, thereby desaturating the torque). We entertain here the question 
of this asymptotic value, and we derive torque formulae which give the total torque as a function 
of the disk's viscosity and thermal diffusivity. The horseshoe drag features two components: one 
which scales with the vortensity gradient, and one which scales with the entropy gradient, and which 
constitutes the most promising candidate for halting inward type I migration. Our analysis, which is 
complemented by numerical simulations, recovers characteristics already noted by numericists, namely 
that the viscous timescale across the horseshoe region must be shorter than the libration time in order 
to avoid saturation, and that, provided this condition is satisfied, the entropy related part of the 
horseshoe drag remains large if the thermal timescale is shorter than the libration time. Side results 
include a study of the Lindblad torque as a function of thermal diffusivity, and a contribution to the 
corotation torque arising from vortensity viscously created at the contact discontinuities that appear 
at the horseshoe separatrices. For the convenience of the reader mostly interested in the torque 
formulae, section 8 is self-contained. 

Subject headings: Planetary systems: formation — planetary systems: protoplanetary disks — Accre- 
tion, accretion disks — Methods: numerical — Hydrodynamics 



1. INTRODUCTION 

Ever since the discovery of the first extrasolar plan- 
ets, many efforts have been undertaken to account for 
the statistics of their orbital properties, and to link 
these statistics to the properties of the protoplanetary 
disk. In particular, the so-called planetary population 
synthesis models correspond to a class of essentially 
one dimensional models, in which the physics of the 
disk is treated so as to include, in a simplified man- 
ner, as many relevant physical processes as necessary, 
and in which the growth and migration of randomly 
sorted pl anetary embryos is tr a cked until the d i sk dis- 
appears (IThommes et al.l 120081: IMordasini et al.l [2 009a: 
Ida fc Lin! l2008bllal. I2004bllal. 120051: lAlibert et all 1200a 
Thommes et all 120071 : IThommes fc Murray! 12001) . The 
statistical properties inferred from these models can be 
confronted to observational data, and thus allow to con- 
strain the underlying disk physical properties or the 
physics of the interaction of the embryos with the disk. 
One drawback of these models, however, is that they 
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have to rely on tidal torque expressions, since the pro- 
cesses that give rise to exchange of angular momentum 
between a planetary embryo and the disk cannot be cap- 
tured by a one dimensional analysis. 

A problematic phase arises in these models, that corre- 
sponds to the migration of low mass objects (the type I 
migration phase). Stand ard torque es tim ates, based on 
the lin ear expressions of IWarl (fT997h or iTanaka et all 
(2002), yield a rapid flush of the low mass cores to the 
innermost parts of the protoplanetary disk. Several au- 
thors have circumvented this problem by arbitrarily low- 
ering the type I migration estimate by an ad hoc factor / 
typically in the range / ~ 0.01 — 0.1, so as to allow a rea- 
sonable fraction of the initial cores to remain as sizable 
distances of the star over the disk lifetime (llda fc Lin! 
I2008al IMordasini et alJ[2009H ). 

The theory of the tidal interaction of low mass plan- 
etary o bjects with the disk has recently ch anged, how- 
ever, as lPaardekooper fc Papaloizoul (|2009aT ) have shown 
that one of the tidal torque components, the so-called 
corotation torque, eventually becomes non-linear at all 
planetary masses. This breakthrough re nders obsolete 
the sy stematic use of the linear estimate of ITanaka et al.l 
((2002T) for low mass planets, and the torque expression 
that should be use d is the sum of the linear Lindblad 
torque expression of lTanaka et al.l (|2002l ). and of the non- 
linear corotation torque. It corresponds to the torque 
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exerted by the coorbital region on to the planet. Only in 
the limit of low mass and large diffusion, does the corota- 
tion torque approaches its linear estimate. Far from this 
limit, the non- linear corotation amounts to the horse- 
shoe drag, which corresponds to the torque exerted by 
fluid elements undergoing a horseshoe like libration in 
the coorbital region, in the frame corotating with the 
planet. 

The horseshoe drag bears a number of similarities with 
the linear corotation torque: its scaling with all parame- 
ters (planet mass, disk surface density, disk aspect ratio, 
etc.) is the same, but it turns out to have in general a 
larger value. Another important difference is its long- 
term behavior: the horseshoe dynamics tends to redis- 
tribute angular momentum between the horseshoe region 
and the planet, but the total angular momentum content 
of the system composed of the planet and its coorbital 
region remains constant in an inviscid disk, if one dis- 
regards what is carried away by the wake. As a con- 
sequence, the horseshoe drag tends to zero on the long 
term, and the total angular momentum ever exchanged 
between the disk and a planet initially "switched on" in 
the unperturbed disk is simply the difference of angular 
momentum content of the horseshoe region between its 
final state and its initial, unperturbed state. This pro- 
cess, known as saturation, occurs on a relatively small 
time scale, typically of O(10 2 ) orbital periods for Earth- 
sized objects. It is in particular much shorter than the 
disk lifetime and the migration time scale itself, hence 
models that follow the migration of a core up to the 
disk dispersal should consider the long term value of the 
horseshoe drag (which we also call the asymptotic or sat- 
urated value) , rather than the value that it displays soon 
after the introduction of the planet in an unperturbed 
disk (which we call the unsaturated value). As stressed 
above, if there is no process that allows the exchange 
of angular momentum between the horseshoe region and 
the rest of the disk, the asymptotic horseshoe drag is 
simply zero. If the disk has a non- vanishing viscosity, 
however, the viscous friction at the boundaries of the 
horseshoe region (the separatrices) enables a steady flux 
of angular momentum across the horseshoe region, and 
the asymptotic value of the horseshoe drag may remain 
finite. It is the purpose of this paper to provide the 
asymptotic value of the torque. 

An estimate of the sat urated horses hoe drag already 
exists in the literature (jMassetl 120011 ). but it applies 
exclusively to isothermal (or barotropic) disks, that is to 
say to disks in which the pressure is solely a function of 
the density. It has recently been shown that the horse- 
shoe dynamics is more complex in non-isothermal flows, 
which obey an energy equation rather than a ba ro tropic 
closure relation (IPaardekooper fc Mellemal 120061: 
Baruteau fc Massetl 120081: IPaardekooper fe Papaloizoul 
2008t iMasset fe Casolill2009t IPaardekooper et al.l l2009) . 

In any case, the key quantity that determines the value 
of the horseshoe drag is the vortensity, or potential vor- 
ticity. It is the value of the vertical component of the 
vorticity divided by the surface density of the disk. Re- 
gardless of the equation of state, the horseshoe drag al- 
ways reduces to a simple integral involving the vorten- 
sity distribution across the horseshoe region. In a non- 
isothermal flow, baroclinic effects act as source or sink 
terms for the vortensity, which in turn has a strong im- 



pact on the horseshoe drag. These effects have been 
found to be localized exclusively at the sep aratrices of 
the horseshoe region (|Masset fc Casolil 12009) . and they 
scale with the entropy gradient in the disk. As entropy 
follows the horseshoe streamlines and tends to be flat- 
tened over the horseshoe region by phase mixing, the en- 
tropy gradient that persists over long time scales across 
the horseshoe region depends on diffusive processes that 
affect the evolution of entropy (such as thermal diffusion, 
and in a more general manner all the processes that act as 
sink or source terms in the energy equation). The asymp- 
totic torque expression therefore involves two main dis- 
sipative processes: viscosity, and thermal diffusion, and 
it depends on two gradients: the gradient of vortensity, 
and the gradient of entropy. 

As a point of terminology, we mention that the 
component of the horseshoe drag which scales with 
the vortensity g radient has been c alled vortensity- 
relate d torque (IPaardekooper et all 120091) . or bulk 
term ()Masset fc Casolil I2009D . whereas its additional 
non-barotropic component has been cal led entropy- 
related torque (| P a ardckoo per et al.l 120091 ) . edge term 
( Mas set fc Casolil 12009ft or adiabatic torque excess 
(|Masset fc Casolil 120091 iBaruteau fc Massetl 12008) . We 
will use indifferently these denominations, but we shall 
reserve the last one for the unsaturated horseshoe drag, 
for which the dissipative processes are unimportant and 
the flow can be considered as adiabatic. 

This paper is organized as follows: section [2] presents 
the framework and introduces the notation and governing 
equations, section [3] presents a reduced model of the flow 
that retains its essential characteristics for the horseshoe 
dynamics, section [4] presents a derivation of the horse- 
shoe drag expression that we shall use to get the asymp- 
totic estimates. This expression is strictl y equivalent to 
the expression of IMasset fc CasoU ((2009) , in which the 
edge term is now incorporated into the integral over the 
vortensity, which is singular at the separatrix. In sec- 
tion [5j we present the set of full hydrodynamical simu- 
lations that we have performed with the FARGO code, 
and we compare their results with those of the reduced 
model of section G3 In section [6j we work out suitable 
relationships for the steady state flow, and infer the cor- 
responding asymptotic horseshoe drag. Our analysis is 
exact in the limit of small viscosity or small thermal dif- 
fusion, but its overall accuracy is found to be satisfac- 
tory up to the values for which the torque remains fully 
unsaturated. For the convenience of the reader mainly 
interested in the results and their implementation, sec- 
tion [5] is self-contained (the notation required to follow 
that section is given in Tab. [IJ. It provides the torque 
value as a function of viscosity and thermal diffusivity 
for low-mass planets subject to type I migration. Sec- 
tion |9] provides additional discussion, and we draw our 
conclusions in section [TU1 

2. PREREQUISITE 

2.1. Framework 

A fundamental premise of our analysis is that a tur- 
bulent protoplanetary disk can be adequately modeled 
by a laminar, non-magnetized disk, in which diffu- 
sion (of energy and angular momentum) is described 
by phenomenological diffusion coefficients. Turbulence, 
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whi ch likely results from th e magnetorotational instabil- 
ity (|Balbus fc Hawlevl [l991'). acts at exchanging angular 
momentum across the horseshoe separatrices. We as- 
sume that it does so in a diffusive manner on the radial 
scale of the horseshoe region. It is not clear whether this 
is the case for the low or intermediate mass planets con- 
sidered here, for which the horseshoe zone width is at 
most of the order of the disk thickness, that is to say 
of same order than the largest scale of the turbulence. 
Since this is a totally unexplored regime, which awaits 
dedicated three dimensional MHD calculations, we shall 
restricts ourselves to this assumption so that it can be 
modeled by resorting to the Navier Stokes equation with 
an effective kinematic viscosity. Little is known on the 
validity of this assumption in the context of horseshoe 
drag saturati on. Nevertheless, we me ntion the encourag- 
ing results of iBaruteau fc L in (2010), who have consid- 
ered an inviscid two-dimensional disk subject to stochas- 
tic forcing, and found saturation properties of the horse- 
shoe drag very similar to those in a laminar disk. 

In this picture the details of the properties of turbu- 
lence should not matter. The scale at which the turbu- 
lent cascade dissipates (resistive or viscous) may impact 
the saturated state of th e turbulence, i.e. the value of 
the effective viscosity (see lFromang et aTll2010l and refs. 
therein). However, provided the correct value of the ef- 
fective viscosity is used, the time averaged distribution 
of vortensity within the horseshoe region should be in- 
sensitive to the statistical properties of the turbulence. 

We note that all what matters to assess the asymp- 
totic horseshoe drag is how much angular momentum 
exchange subsists at large time between the disk and the 
horseshoe region, considered here in a wide meaning as 
everything enclosed with in the separatrices, including the 
planet (e.g. iMa ssct 2001). For this reason, the horseshoe 
drag is always linked to the flow properties at the sepa- 
ratrices. The internal details of the coorbital flow, such 
as the existence of tadpole streamlines or the topology of 
the flow in the vicinity of the planet, do not impact the 
horseshoe drag. Stated differently, since the whole set of 
librating material is trappe d in the coorbital region an d 
migrates along with it (e.g. Masset fc Papaloizoull2003l ). 
the forces internal to the horseshoe region are internal 
to the migrating system and therefore do not affect its 
migration rate. 

Turbulence also yields a stochastic torque that trig- 
gers a random walk of the semi-major axis of the planet 
(|Nelsonl[2005h . which comes in addition to the system- 
atic drift due to the wake. As a consequence the prob- 
ability density of the p lanet's semi-major a xis obeys a 
diffusion-advection law (| Johnson et al.ll2006f) . When av- 
eraging over a sufficient amount of time, the systematic 
drift overcomes the stochastic effects. Although it is still 
unclear whether the drift so averaged amounts to the 
drift obtained from a laminar analysis ([Nelson 2005|), we 
assume here that this is the case, and that the long-term 
effects of turbulence exclusively amount to regulating the 
degree of saturation of the horseshoe region, from which 
we infer laminar torque values which dictate the time 
averaged migration rate. 

In a similar vein, a number of different processes, in ad- 
dition to turbulence itself, lead to a diffusion of tempera- 
ture in the disk, such as thermal diffusion or radiative dif- 
fusion, while cooling through the disk's photosphere and 



heating by turbulence also contribute to the energy bud- 
get. For the sake of simplicity, we shall assume that the 
time evolution of the fluctuations of temperature with 
respect to the unperturbed state can be described by a 
unique diffusion coefficient, to which we refer as the ther- 
mal diffusivity, which is an effective diffusion coefficient 
much as the effecti ve viscosity. This is a rather stan - 
dard practice (e.g. iPaardekooper fc Papaloizoul 120081) . 
We shall see that all what matters to determine the de- 
gree of saturation of the entropy related torque is the 
entropy difference on both separatrices upstream of the 
U-turns, and that the only parameter that determines 
this difference is the ratio of the thermal timescale to 
the horseshoe libration timescale. Our simplifying as- 
sumption therefore retains all the physics relevant to the 
description of the saturation of the non-isothermal horse- 
shoe drag. 

Finally, most of our study (analytical and numerical) is 
performed assuming a two-dimensional disk, for reasons 
of tractability and computational cost. In section 17.2.11 
we will try and generalize to three dimensions in the 
framework of stacked two dimensional layers. 

2.2. Notation 

In addition to the notation defined below, the nota- 
tion used in section [5] is presented in table [TJ for the 
convenience of the reader mainly interested in the re- 
sults presented in that section. We consider a planet of 
mass M p orbiting a star of mass Af*, on a fixed circu- 
lar orbit. This planet is embedded in a gaseous proto- 
planetary disk, such that the planetary orbit is coplanar 
with the disk and prograde. We use P to denote the 
vertically integrated pressure, and s to denote the gas 
entropy, which reads: 



P 

£7' 



(1) 



where E and 7 are defined in Tab. [T] 

We identify a location in the disk by its distance r 
to the star and its azimuth cf> with respect to the planet. 
The gas has a radial velocity v r and an azimuthal velocity 
U0 in the frame corotating with the planet (hence its 
angular frequency in an inertial frame is Q — v^/r + 
flp). The radius of corotation r c is the location in the 
unperturbed disk where the material has same angular 
frequency as the planet: ft(r c ) = Q p . We will make use 
of the distance x to corotation: x = r — r c . We also 
consider the disk pressure scaleheight H = c s /Q = rh. 

We denote with a subscript the value of a variable 
in the unperturbed disk, and with a subscript 1 the per- 
turbation of this variable. Whenever we refer to the un- 
perturbed value of a variable at corotation, we use a c 
subscript. 

Using the notation of IBaruteau k, Massed ()2008[ ). we 
define the gradient of the inverse of vortensity across 
corotation as: 



V = 



dlog(S /wo) 



dlogr 



r c gK 
£ r dr 



(2) 



where a is the index of the surface density power law 
(see Tab. [T|), and where u> = (l/r)9 r (r 2 J7) is the verti- 
cal component of the flow vorticity. We will also make 
use of the first Oort's constant A — (l/2)rdfl/dr, which 



4 



F. S. Masset & J. Casoli 



Table 1 

Main notation used in this work. 



Name 


Notation 


Comment 


Planet orbital frequency 


^ p 










Planet to star mass ratio 


g 


q = Mp/M* 


Disk surface density 


s 


S = S c (r/a)" Q 


Vertically averaged temperature 


T 


T = T c {r/a)-P 


Flaring Index 


/ 


/ = (l-/3)/2 

h QC 


Disk aspect ratio 


h 


h = ho(r/a)f 


Vortensity gradient 


V 


V = 3/2 - a 


Adiabatic index 


7 




Entropy gradient 


S 


7 7 
The entropy scales as r" 10 






Horseshoe half width 


Is 





Note. — The entropy gradient £ defined by Paardckoopcr ct al. (2009) 
is equivalent to yS. 



quantifies the amount of shear in the unperturbed flow, 
and of the second Oort's constant B = (1 / '2r)d(r 2 Vl) / dr , 
which is half the vertical component of the unperturbed 
flow's vorticity. Instead of the vortensity w defined by 
w = w/E, we oftentimes consider its inverse I = 1/w, 
that we call the load for the sake of definiteness: I — E/w. 
We define the gradient of entropy across corotation as: 



S = 



1 d log s 



7 d log r 

We shall make use of the dimensionless variables 
L = (l — l c )/l c and S=(s-s c )/s c . 



(3) 



(4) 



We denote the slope of these variables in the unper- 
turbed disk respectively with L'^ and S'^, where the 
prime stands for the derivative with respect to x. We 
have the following relationships: 



V 

a 



(5) 
(6) 



The reason for the oo symbol in the expressions above 
will appear in section 13.31 They represent the large 
scale gradient in the unperturbed disk, and therefore the 
asymptotic value of the slope of the load or entropy as 
one goes away from the coorbital region. 

2.3. Governing equations 

The governing equations of the flow are the equation 
of continuity, the Navier-Stokes equations and the energy 
equation, together with the closure relationship provided 
by the equation of state, which is that of an ideal gas: 



P = KZT/n, 



(7) 



where 1Z is the ideal gas constant and p, the mean molec- 
ular weight of the gas. The equation of continuity reads, 
in the frame corotating with the planet: 

5 t E + id r (Erw r ) + ^d (Et; ) = 0. (8) 



The Navier-Stokes equations read, respectively in the ra- 
dial and azimuthal directions: 



dtV r +v r d r v r +—d t f > v r —rVL p —2Q, p v t f > - 



and: 



(9) 



(10) 



(11) 



where D t = d t + v r d r + ^-8$ and 
3 = r 2 n 

is the specific angular momentum. In Eqs. ([9]) and (|10p . 
f r and frf, are respectively the radial and azimuthal com- 
ponent of the viscous force per unitary surface. These 
viscous forces are derived from the viscous stress tensor 
as follows: 

j = 1 d{rr rr ) | 1 dr ri p r 00 
r dr r d<p r 

U = 5 1 oT ' ' i L6 ) 

r or r 0(f> r 

where the components of the viscous stress tensor are: 

2 



-2r]D r 



77V. v 



where 



T00 = 2?7£)00 - -??V.v 



dv r 1 dv<p v r 

or r od> r 



dr \ r J r 



1 dv r 



(14) 

(15) 
(16) 

(17) 
(18) 



and rj = E/v is the vertically integrated dynamical vis- 
cosity coefficient. The internal energy density is e = 
p/("f — 1), and the energy equation reads: 



ED, 



f- 

Ve 



-pV.v- V.F, 



(19) 
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where the last term of the right hand side accounts for 
the energy diffusion, and will be given explicitly in sec- 
tion EL51 



then be recast as: 



2.4. Time evolution of the load 

The evolution of the load or vortensity that the ma- 
terial undergoes in the coorbital region corresponds to a 
slow transformation which takes place while the fluid el- 
ements describe essentially circular trajectories, far from 
the planet, between successive horseshoe U-turns. We 
therefore neglect the azimuthal derivatives in the equa- 
tion governing the evolution of the load. Eq. (| 10[) can 
then be recast as: 



d t j + v r rcj = -^<9 r (r 2 T0 r ), 



(20) 



where we have used (fTTj) . (JT5J) and (TT5T) . Deriving with 
respect to r, using Eq. ((§} and (fl"8|). we are left with: 



D t w 



d r (Evr 3 d r {l) 



(21) 



We assume that v is a power law of radius, such that the 
inner bracket of Eq. (|21[) is flat in an unperturbed disk, 
which amounts to assuming that there is no radial drift 
of material. We mention that, in some of our numerical 
simulations, we have precisely adopted such a prescrip- 
tion (see section 15.3.21) . We can thus write, keeping only 
the highest order derivatives of the perturbed quantities: 



D t w ps vaz 2 W\ — 2vw 



(22) 



where Slo and So are respectively the angular velocity 
and surface density of the unperturbed disk, whereas fii 
and Si are the perturbations of azimuthal velocity and 
surface density. In Eq. (|22|) we have made use of the first 
order expansion: 



Wl ~ r s 1? 

and of the relationship, specific to Keplerian disks: 

r Oo = — 3 . 

r 

Recasting Eq. ([22]) in terms of the load, we obtain: 



D t l*>vd%l + 2vh^ 



d 2 ^ 



(23) 



(24) 



(25) 



where we have used the relationships d% 2 h ~ d^ 2 l and 
d^ 2 £i ps d^ 2 S. The typical value of the perturbed load 
in the horseshoe region is indeed 0(l c x s /a), whereas it 
varies radially over a length scale at most equal to x s . 
The same applies to the surface density, except in a 
barotropic disk: we note that if the radial scale of the 
disturbances in the disk is small with respect to the pres- 
sure scale length (i.e. with respect to the disk thickness), 
then, in a barotropic disk, the excitation of evanescent 
pressure waves renders the relative perturbation of sur- 
face density much sm aller than the relative perturbation 
of vortensity or load (jCasoli &: M assct 20(53) , so that the 
perturbation of vortensity can be put solely on the ac- 
count of the perturbation of vorticity, and Eq. (|25|) can 



Dtlmvd^l. 



(26) 



On the contrary if the pressure length scale is shorter 
than the radial scale of the disturbances, then the acous- 
tic spread of the perturbations does not significantly alter 
the perturbations of surface density, so that the relative 
perturbation of load is equal to the relative perturbation 
of surface density. In that case, Eq. (|2~5T) yields: 



D t l « ivdiil. 



(27) 



In the case that we contemplate in this work, that of a 
low-mass planet embedded in a disk, for which the half- 
width of t he horseshoe regio n is smaller than the disk 
thickness (|Masset et all 120061 ). the regime that prevails 
corresponds to Eq. (|26)l . Eq. (|27l) would be valid for very 
massive protoplanets, but it would be of little interest, 
since those clear a gap in their coorbital region, thereby 
shutting off the corotation torque. 

We note that in a non-barotropic disk, in which there 
can be entropy waves and therefore contact disconti- 
nuities, the last term of Eq. (|25|) may no longer be 
negligible. This is in particular the case at the outer 
edge of the horseshoe r egion, downstream of the U-turns 
(jMasset fe Casohll2009D . We shall therefore use Eq. (gSJ 
for barotropic disks, and Eq. (|2"5"|) for disks with an en- 
ergy equation. 

2.5. Time evolution of the entropy 
Using Eq. ([7]) we can transform Eq. ([T9"]) into: 

TZYIT 

— Alogs = -V.F, (28) 

(7- 1)M 

where F is the heat flux, given by: 

F = -fcVT, (29) 

where k is the thermal diffusivity. Assuming that it can 
be regarded as constant over the length scales over which 
T varies, Eq. (|2"51) can be recast as: 



TZY.T 
(7-1)/* 



D t log s = kAT. 



(30) 



We now note that, for the case we are interested in, that 
of a horseshoe region bound by contact discontinuities, 
the variations of temperature occur on a much smaller 
scale than those of pressure, which occur on the disk 
scale height. We therefore make the following simplifying 
assumption, using Eq. ([I]): 



AT _ 1 As 

T 7 s 

Using Eqs. (|50"f and (pITj) . we are led to: 

D t s — kAs, 

where 

7 — 1 kfj, 

The simplifying assumption underlying Eq. (1311) has sev- 
eral important properties: 



(31) 



(32) 



(33) 
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• Our prescription for thermal diffusion has no im- 
pact on the acoustic waves, since the diffusion term 
scales with the perturbations of entropy, which van- 
ish in acoustic waves. It has therefore no impact 
on the differential Lindblad torque, and allows to 
disentangle in a clean manner the variations of the 
corotation torque due to variations of thermal dif- 
fusion. 

• Eq. ([32|l shows that the evolution of entropy does 
not depend on other variables, if we regard E as 
a constant, and if we consider that the streamlines 
do not depend on the perturbation that builds up 
in the coorbital region. Finding a steady state so- 
lution for the coorbital flow can therefore be done 
in two subsequent steps: (i) one first seeks a steady 
state solution for the entropy, which does not re- 
quire the knowledge of any other variable, (ii) one 
then seeks a steady state solution of the vorten- 
sity, as the latter depends on the entropy field, in a 
manner that shall be specified in the next section. 

We comment that we have undertaken a side study of the 
variation of the Lindblad torque as the thermal diffusiv- 
ity varies (see section [77Ll) using the alternate form of the 
diffusion term (which then scales with the Laplacian of 
the temperature, rather than the entropy). The transi- 
tion of the Lindblad torque from the adiabatic limit to 
the isothermal limit (the respective torque values differ 
by a factor 7) occurs for values of the thermal diffusion 
much larger than those relevant for the saturation of the 
entropy related corotation torque. 

3. A SIMPLIFIED MODEL OF THE FLOW 

3.1. Simplifying assumptions 

Owing to the complexity of the coorbital flow of an em- 
bedded planet, we have to make a number of simplifying 
assumptions which render tractable the task of finding 
the expression of the steady state horseshoe drag. We 
anticipate that these assumptions do not impact signifi- 
cantly the horseshoe drag expression. We will check this 
in section 15^41 

• We separate the time scales of the horseshoe U- 
turns and that of the drift between two successive 
U-turns (i.e. half the libration time scale). The 
former is at leas t one order of magnitude shorter 
than the latter (jBaruteau fe Massed l2008h . This 
allows us to consider that the horseshoe U-turns 
are performed instantaneously. 

• Similarly, we consider that the azimuthal interval 
over which the U-turns take place is small com- 
pared to 27r. This amounts to considering that the 
U-turns are squeezed against the axis joining the 
star and the planet. 

• We assume that the vortensity (or load) is con- 
served during the U-turns, except for the stagna- 
tion streamline (separatrix) in the non barotropic 
case, at w hich a singular amount of vortensity 
is created (|Masset fe Casolil l2009f ). Note that 
some vortensity is cr eated during the U-tur ns in a 
non barotropic flow (|Masset fc Casolill2009() for all 




Figure 1. Sketch of the flow and of the load profile in front and 
behind the planet, within our simplified framework. The graphs 
on the left and right side present a sketch of the front and rear 
load profiles, respectively. The thick arrows in these graphs rep- 
resent the singular load on the downstream separatrix. The light 
shaded area enclosed within the horizontal dashed lines shows the 
horseshoe region. The functions L F (x) and L R (x) are even for 
x E [— x s ,x s ]. The central set of arrows depicts the Keplerian 
flow, for which v y = 2Ax. 

streamlines, but this has no impact on the horse- 
shoe drag, hence we disregard it in this analysis. 
Also, the vortensity is ill-defined at the downstream 
separatrix, where there is both a contact discon- 
tinuity and a vorticity sheet. Nevertheless, the 
contact discontinuity has a s urface density jump 
which is first order in x s /a (|Baruteau fc Massed 
2008), hence to lowest order mx B /a the singularity 
of vortensity can be defined as Av < j,S(x — x s )/T,o, 
where Av^ is the jump of azimuthal velocity across 
the downstream separatrix. 

• We assume that the streamlines can be regarded 
as circular between two U-turns. This amounts to 
neglecting the radial velocity in the coorbital re- 
gion, and to disregarding the behavior induced by 
the indirect term of the potential. 

• We assume that the U-turns are symmetric, i.e. 
that a streamline which has a radius r prior to a U- 
turn is mapped to a streamline with radius 2r c — r, 
its distance to corotation remaining therefore \r — 

To\. 

A consequence of this set of assumptions is that the path 
followed by a horseshoe fluid element in the (y, x) plane 
is rectangular (where y = rcf>), and that the load and en- 
tropy are even functions of x over the interval [—x s , x s ] 
in <p = (front U-turns) and in (f> = 2n (rear U-turns), 
except for the possible creation of a singular amount of 
load in x — ±x s . Hereafter we will write with a super- 
script F (R) the quantities in y = + (y = 2ira~) so as to 
remind that they are evaluated in front (at the rear) of 
the planet. A generic situation for the load is represented 
in Fig. [IJ 

3.2. Governing equations of the reduced model 

In the framework considered in the previous section, we 
consider a restricted model of the flow, in which we have 
only one scalar field in the barotropic case (the load) or 
two scalar fields in the case with an energy equation: the 
load and the entropy. In the barotropic case, the variable 
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L obeys an equation derived from Eq. (|26[) . which reads 
in the framework of the reduced model: 



d t L + 2Axd y L 



(34) 



In the case with an energy equation, the variable L obeys 
an equation similar to Eq. (l25j) . The radial scale for 
the variation of pressure is much larger than the scale of 
variation of entropy or surface density, so we are led to: 



dtL- 



2Axd v L = vL" - —S". 

7 



(35) 



In this case, we also need to consider the time evolution 
of the entropy. Eq. (|3"2")l yields: 



d t S + 2Axd v S = kS". 



(36) 



Eqs. (|3"4"1) to (|3"6l amount to considering exclusively the 
radial diffusion of these quantities, and to assuming that 
the underlying flow is the unperturbed Keplerian flow, 
rectified on a Cartesian slab. 

3.3. Boundary conditions of the reduced model 

The boundary conditions of this reduced model are pe- 
riodic in y outside of the horseshoe region, with a period 
equal to 2na, that is to say: 

L F (x) = L R (x) and S F (x) = S R (x) for any x such that \x 

(37) 

whereas within the horseshoe region the following rela- 
tionships hold: 

L F {x) = L F {-x) +L 6(x + x s ) for any x € [-a; s ,$8) 
L R (x) = L R (-x) - L Q S(x - x s ) for any x £ [0, x s ] (39) 
and 

S F (x) = S F {-x) for any x G [-x s ,x s ] (40) 
S R (x) = S R (-x) for any x £ [-x s , x s ] (41) 

Eqs. (|3"5)) and (131?)) are meant to account for the horseshoe 
U-turns, in which the load is conserved, except at the 
downstream separatrices where a singul ar amount ±Aq = 
LqI c i s created, as has been shown by iMasset fc Casoli 
(2009), and as we shall discuss more in detail in section |4j 
The entropy obeys similar equations, except that it is 
conserved during the U-turns, even at the separatrix. 

In the x direction, boundary conditions are set for 
| x | — ^ oo such that the slopes of the fields L and S are 
that of the unperturbed flow: 

lim d x L(x,y) — L'^ for any y £ [0,2™] (42) 

x— >±oo 

and 

lim d x S(x, y) = S'^ for any y £ [0, 2na] (43) 

x— ^±oo 

3.4. Symmetry of the reduced model 
The reduced model exhibits the following symmetry: 

L(x,y) = —L(-x,2na — y). (44) 

while a similar property holds for the entropy field. 
The system is set with this symmetry at t = 0, since 
L(x,y) oc x in the unperturbed flow, and the governing 



equation (|34l) as well as the boundary conditions speci- 
fied by Eqs. (|37l) to (|39j) ensure that Eq. gU) holds at 
any time. As a consequence we have: 



and 



L (—x) = —L F (x) for any x 
S R (-x) = -S F (x) for any x 



(45) 
(46) 



3.5. General considerations on the steady state of the 
reduced model 

The scalar field S is not subject to any feed back from 
L: it evolves by itself. On the contrary, the singular 
vortensity ±Ao created at the downstream separatrices 
is linked to a difference of the value of t he entropy at the 
upstr eam front and rear separatrices ( Masset fc Casoli 
2009). We have, using Eq. (62) of IMasset fc Casoli 
([20091) : 

_ SqAi; _ I\ 
(2Bf 16a\A\x 2 s B 2 ' 1 ' 

where T\, the adiabatic torque excess, has been worked 
out bv IMasset fc Casolil (|2009fl [see e.g. their Eq. (96)], 
and where Av is the azimuthal velocity jump across the 
downstream separatrix. In an unsaturated situation, Fi 
(and therefore Ao) depend on the unperturbed entropy 
gradient S in the disk, which imposes the entropy diffcr- 
>caee at the upstream separatrices. In a saturated situa- 
tion, the latter must be substituted by the effective en- 
tropy gradient, which depends on the actual value of the 
entropy difference at the upstream separatrices (rather 
than the unperturbed one): 



S' = - 



1 a s F (x s ) — s (— x s ) 1 a 



7 «c 



2lC a 



--S F (x s ), (48) 
7 x s 



where we have used the symmetry property of Eq. (|46[) . 
The steady state situation is therefore determined as fol- 
lows: 

• A steady state solution is sought for the entropy 
field 5. 



• Using Eq. (|48|) . this solution is used to evaluate Ti 
and therefore A , the production of vortensity at 
the downstream separatrix. 

• Once Ao is known, one can determine a steady state 
solution for the load field L. 

• The latter can then be used to evaluate the horse- 
shoe drag as explained in section |4j 

3.6. Implementation of the reduced model 

The reduced model presented above can be easily im- 
plemented as a simple grid-based code. It allows to get 
values of the asymptotic horseshoe drag at low com- 
putational cost. Namely, we solve the evolution equa- 
tions f34|) and (|36|) on a mesh subdivided in iV x x AL 
zones. The advection of the field (either load or entropy) 
can be ensured by any upstream TVD scheme (in our 
case we use a monotonized centered slope limiter). The 
timestep is limited by the Courant criterion, which reads 
here: 



At = Cq min[(At7 2 + At^ 2 + At^y 1 ' 2 ] 



(49) 



F. S. Masset & J. Casoli 



LOW VISCOSITY 



bv lMasset fc~Caso li (2009) can be recast as: 




Figure 2. Still frame of the animation, that shows the load dis- 
tribution at the end of the calculation (bottom) and the horseshoe 
drag as a function of time, in arbitrary units (top). The viscosity 
is low in this case, so that the load distribution eventually becomes 
fiat. The animations shows two additional cases: a high viscosity 
case, in which the torque levels off to a non-vanishing value, and 
a case with an energy equation (i.e. the appearance of singular 
load sheets at the downstream separatriccs) in which the torque 
oscillates and eventually decays to a rather low value. 



where Aii = Ax 2 /(2v), At 2 = Ax 2 /(2n), At 3 = 
Ay/\2Ax\, Ax and Ay being the mesh resolution respec- 
tively in x and y (Ax <C Ay). In Eq. (|49|) we chose 
the Courant number Co = 0.9, and the min operator is 
meant to be taken for all the zones. The quadratic sum 
that features in Eq. (l49j) is similar to that in Eq. (74) of 
iStone fe Normanl (|1992f ). 

The diffusion either of load or entropy is implemented 
by a second order explicit operator, in the x direction 
only. The corresponding substep therefore reads, for the 
load: 



jb _ TO. 



vAt- 



2L ?i 



(50) 



where i and j are the zone indexes respectively in x and 
y, and where the a and b superscripts denote the load 
respectively before and after the diffusion substep. 

The mesh is designed so that the corotation and both 
separatrices fall on the interface between adjacent rows 
of zones. The horseshoe region therefore covers an in- 
teger and even number of rows. One column of ghost 
zones is used on each side in y. For the zones of the in- 
ner or outer disk, standard periodic boundary conditions 
are applied in order to set the values in the ghost zone. 
For the rows which are within the horseshoe region, the 
boundary conditions differ, so as to mimic the U-turns. 
An animation of the time dependent load distribution is 
available in a separate file. Fig. [2] shows a still frame ex- 
tracted from this animation, which also displays the cor- 
responding horseshoe drag value. The manner in which 
the latter is estimated is detailed in the next section. 

4. HORSESHOE DRAG EXPRESSION 

We use the horses hoe drag expression given by 
iMasset fe Casolil ([2009I ). which we transform so that the 
entropy related torque (corresponding to the last term of 
this equation) feature the load produced during the U- 
turns. The horseshoe drag expression given at Eq. (37) 



r = r F + r R , 



(51) 



where the superscripts F and R respectively stand for 
the front and rear contribution to the horseshoe drag, 
which read respectively: 



T F = f *l + Aj dG-^Af AG 



1-AjodG - §A Jo s AG, 



(52) 
(53) 



where G is defined at Eq. (16) of lMasset fc Casolil (|2009f ). 
and where Ajo = AaBx is the jump of specific angular 
momentum of fluid elements that execute a U-turn from 
r = a + x to r = a — x, while Aj'q refers to the same quan- 
tity for x — x s . Note that we have split the torque excess 
(trailing term of Eq. (37) in MC09) into a front part and 
a rear part, which are equal for symmetry reasons. Using 
Eqs. (19), (37) and (62) of lMasset fc Casolil (pool ), one 
can recast the front contribution to the horseshoe drag 
as an integral over < x < x s , corresponding to the 
streamlines which undergo a horseshoe U-turn towards 
smaller radii in front of the planet: 



4\A\Bl F (x)Aj .xdx + 2| A\B A A j%x s , (54) 



r F = 



By convention, Eq. (J54J) and all the torque expressions 
that we shall write correspond to the torque exerted by 
the disk onto the planet. From the definition of Ajo 
it follows that it is a positive quantity on the interval 
considered, and that T F is also positive, as expected since 
it stems from material in front of the planet. In Eq. (1541) . 
the trailing term come s from the torque e xcess expression 
given by Eq. (37) of IMasset fc Casolil ([2009D . that we 
have transformed using Eq. (|4"T)) . We recast Eq. ([54]) as: 

r F = 16\A\B 2 a [ l F {x)x 2 dx + 8\A\B 2 aA Q x 2 s , (55) 



We transform Eq. ([53)) using Eq. ([3"5]) , which allows to 
embed the singular contribution within the integral: 



8\A\B 2 a / l F {x)x 2 dx. 



(56) 



In this equation the integrand is singular in x = —x s , 
that is to say at the edge of the integration domain. It 
is therefore necessary to specify explicitly how to deal 
with this singular contribution. Comparison of Eqs. (|55p 
and (|56p shows that one must take into account the whole 
singular distribution in the integral (as if it was approx- 
imated by a narrow function whose support was entirely 
contained within the interval of integration). Hereafter, 
whenever we shall consider a horseshoe drag-like integral, 
with an integrand singular on an edge (i.e. at a separa- 
trix), we shall implicitly consider that all the singularity 
must be taken into account in the integral. 

The rear contribution to the horseshoe drag can be 
worked out in a similar manner and reads: 



T 1 



8\A\B 2 a / l R {x)x 2 d Xl 



(57) 
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where, again, l R (x), the load at the rear of the planet, 
may contain a singular contribution in x = x s , arising 
from the production of load during the U-turn. The 
whole horseshoe drag therefore reads: 



r = 4\A\BaZ Q [L F (x) - L R (x)]x 2 dx, 



(58) 



This expression can be further simplified, using Eq. (j?5 
as: 

Ft \ 2 



(59) 



\A\BY> a / L (x)x dx, 



Simple particular cases can be recovered from Eq. (|b 

• In the unsaturated barotropic case with a vorten- 
sity gradient V = aL'^, we have L F (x) = Vx/a, 



from which we infer: 

T = 4|A|BS Vx^, 

or: 

r = r n v, 



(60) 
(61) 

which is t he standard barotropic ho rseshoe drag 
expression (|Wardlll991t IMassetll2001h , where 



T = -£ofi 2 a£, (62) 
if we specialize to the case of a Keplerian disk. 

• In an inviscid barotropic disk in steady state in 
the corotating frame, the load is conserved along a 
streamline, hence L F (x) = L (x). Since L R (x) = 
L R {-x) = -L F (x), we infer that L F {x) = for 
any x £ [0, xg]- The horseshoe drag therefore can- 
cels out. It is fully saturated, and the load is profile 
is fiat across the horseshoe region. 

• In an adiabatic flow without vortensity gradient, 
the only contribution to the integrand stems from 
the singularity Ao at the downstream separatrix, 
which leads to: 



?2 ax 2 a A , 



T = 16|A|^ 
which, using Eq. (j47j), yields r — I\. 

5. NUMERICAL SIMULATIONS 



(63) 



In addition to the calculations performed with the re- 
duced model presented in section [3J we have performed 
three series of full hydrodynamical simulations. We de- 
scribe hereafter the code and the different setups that we 
have adopted. 

5.1. Code and setup 

We have u sed the hydroc ode FARGC0 in its isother- 
mal version (iMassetl l2000allbh and in its adiaba tic ver- 
sion (|Bariteau^Masse3l200l Hiintii^ HOOl). The 
FARGO code is a staggered mesh hydro-code on a po- 
lar grid, with upwind tr ansport and a harmonic, sec- 
ond order slope limiter (|van Leerl Il977f ) . It solves the 
Navier-Stokes and continuity equations for a Keplerian 
disk subject to the gravity of the central object and that 
of embedded protoplanets, and the energy equation in 

3 http://fargo.in2p3.fr 



the case of the adiabatic version. It uses a change of 
rotating frame on each ring that enables one to increase 
significantly the time step. The hydrodynamical solver 
of FA RGO resembles the wid ely known one of the ZEUS 
code (jStone fc Norman 1992), except for the handling of 
momenta advection. The Coriolis force is t reated so a s 
to enforce angular momentum conservation (Klcy 1998). 
The mesh is centered on the primary. It is therefore 
non-inertial. The frame acceleration is incorporated in 
the potential indirect term. 

5.1.1. Dissipation properties 

Viscosity — The full viscous stress tensor in cylindri- 
cal coordinates of the Navier-Stokes equations, given by 
Eqs. flU} to ([18]>. is implemented in FARGO. As we 
shall see by comparing our results to that of the reduced 
model, which features only the diffusion of load due to 
the radial shear in Eq. the saturation properties rely 
essentially on the radial shear, and the other components 
of the viscous stress tensor appear to be unimportant. 



Diffusivity — As discussed in section l275l we assume that, 
for our purpose, the evolution of the entropy can be de- 
scribed by a diffusion equation. Since the entropy is not 
one of the primitive variables handled by the hydrody- 
namical solver, the thermal diffusion substep updates the 
internal energy. This substep amounts to solving the fol- 
lowing partial differential equation: 



d t e = e—d r (rd r logs). 



(64) 



No diffusion therefore occurs in the unperturbed disk, 
where s is a power law of the radius. 

5.1.2. Units and common parameters 

A number of parameters are common to all the runs 
that have been performed. These are: 

• The planet to star mass ratio is q = 1(U 5 (this 
translates into a planetary mass of M p = 3.3 M® 
if the central star has a solar mass.) 

• The disk aspect ratio at the location of the planet 
is h = 0.05. 

• As our calculations are two dimensional, we have 
to soften the planetary gravitational potential $ p 
over a length scale e: 



= — 



GM n 



vV 2 + e 2 



(65) 



We performed our calculations with e — 0.5H. This 
value of e is meant to give results representative of 
the dynamics in a slice of the disk at \z\ ~ 0.5H. 
We note that our results depend indirectly on the 
value of e, since, as we shall see, the degree of satu- 
ration of the horseshoe region depends on its half- 
width x s , which is itself a function of the softening 
length. A systematic study of the dependence of x s 
on the softening length, as well as a study of the 
correspon ding unsaturated torqu e, has been per- 
formed bv lMasset fc Casoiil (|2009D . In section OH 
we add up the contributions from slices at all z to 
get a tentative estimate of the three dimensional 
torque. 
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• We strictly remain in the framework either of the 
barotropic a pproximation, or i n the framework of 
the study of IMasset fc Casolil (|2009l ). who assume 
a flat temperature. In both cases, the initial tem- 
perature profile must be flat. In the FARGO code, 
the temperature profile is imposed through the pa- 
rameter / (called flaring index), defined by: 



h = h [ — 

Jo 



(66) 



In a disk with a flat temperature profile, the pa- 
rameter / must be set to 1/2. 

• In the runs with an energy equation, the adiabatic 
index 7 of the gas is set to 1.4. 

5.2. Sets of calculations 

A general situation involves a disk with arbitrary 
vortensity and entropy gradients V g and S g , and arbi- 
trary values of viscosity and diffusivity. In this situation, 
the steady state solution for the entropy, which evolves 
independently of the load, yields the amount of singu- 
lar load L g that is produced at the downstream sep- 
aratrices. The general situation therefore amounts to 
finding a steady state solution L for the load, that ful- 
fills Eq. (|35|) and the boundary conditions specified at 
Eqs. (07]) to (£11 and (gp with L Q = Lg and L'^ = V g /a. 
Owing to the linearity in L of Eq. (|35l) . we notice that 
if we have two solutions of this equation L\ and L 2 , 
which respectively obey the boundary conditions with 
L'^ = Vg/a and A = for L\ (hence S — 0), and 
L 1 ^ = and A = L g for L 2 , then L\ + L 2 is the solution 
sought for the general problem. The field L\ corresponds 
to a barotropic situation (no load is created) with an arbi- 
trary vortensity gradient, while the field L 2 corresponds 
to a non-barotropic situation with a vanishing vortensity 
gradient. 

This allows us to split the study in two distinct parts: 

• A study of the saturation of the vortensity related 
torque, which we perform for a globally isothermal 
disk. This study only requires varying the viscosity. 

• A study of the saturation of the entropy related 
torque, in a disk without vortensity gradient. The 
asymptotic torque value depends on two parame- 
ters: the viscosity and the thermal diffusivity. 

The torque in a general situation can be deduced by 
adding the asymptotic values of the torque in these more 
restricted situat i ons. We nevertheless comment that 
IMasset fc CasoU (pOOl ) found a significantly larger width 
(by up to 10 — 15 %) of the horseshoe region at large 
value of so that the entropy related torque cannot 
solely be put on the account of the production of load at 
the downstream separatrices. This additional complexity 
may render slightly inaccurate the procedure of adding 
up the results of the separate cases. 

From the above considerations, we performed three sc- 
ries of hydrodynamical calculations: one in a globally 
isothermal disk, in which we vary the viscosity, and two 
in a disk with an energy equation, in which we respec- 
tively vary the viscosity and the diffusivity. We sum up 
in Tab. [2] the specific parameters of these series. The 



mesh has a radial extent that ranges from i? m i n to -R ma x, 
evenly spaced in 7V ra( j annuli. Azimuthally, it is divided 
evenly in N scc sectors. The calculations were ran for up 
to either 500 or 10 3 orbits, although those which had evi- 
dently reached their asymptotic torque value before that 
date were stopped manually. 

We note from Tab. [5] that the isothermal calculations 
are performed with a flat surface density profile (hence 
V = +3/2). We do not have any freedom in the choice 
of the disk parameters for the runs with an energy equa- 
tion: since those must have V = and since they must 
have a flat temperature profile, we find that they have an 
entropy gradient S ~ 0.43. Since this value is positive, 
we expect a negative entropy related torque. 

We also note that the small value of the diffusion co- 
efficients adopted in some of our runs raises the possi- 
bility that numerical diffusion due to grid effects could 
be important. Although we have not undertaken a sys- 
tematic study of those, we are confident that they are 
at most of the order of the smallest diffusivity consid- 
ered here (~ 10~ 9 a 2 n p ). Should that not be the case, 
numerical diffusion would wash out the torque depen- 
dence on viscosity or thermal diffusivity, when either of 
these quantities is small. As we shall see below, we get 
a definite, measurable dependence of the torque on the 
diffusion up the smallest value of the diffusion coefficient 
(except when the flow becomes chaotic, but even though 
one can perform time averages and recover a dependence 
on the diffusion). This clearly indicates that the numer- 
ical diffusion present in our runs is negligible. 

5.3. Results 

5.3.1. Isothermal runs 

The results of set 1 are displayed in Fig. [3] An ad- 
ditional inviscid run is shown, whic h disp l ays th e char- 
acteristic sawtooth shape found by Ward (20Q7J) in the 
framework of the reduced model. Higher viscosity runs 
converge much before t = 500 orbits, while those with 
viscosity v < 2.10~ 8 a 2 f2 p still display oscillatory behav- 
ior at that date. Evaluating the horseshoe drag from 
these runs amounts to subtracting the differential Lind- 
blad torque from the total torque. We have noticed that 
the successive minima and maxima of the horseshoe drag, 
at very low but non- vanishing viscosity, are in geometric 
progression, with a common ratio r G] — 1, 0[. Differently 
said, if we denote with 71 and 73 the value of two succes- 
sive minima (maxima) of the total torque and with 72 the 
value of the maximum (minimum) in between, then the 
value of the differential Lindblad torque that we 

infer assuming that (7$ — ^^lr)^ {1,2, 3} are m geometric 
progression, which is: 



AT LR = 



7173 ~ 72 

7i + 73 - 272 ' 



(67) 



turns out to be remarkably independent of the set of suc- 
cessive extrema that we consider. This allows us to give a 
relatively accurate estimate of the differential Lindblad 
torque, as well as an estimate of the asymptotic value 
of the torque at low viscosity (which would otherwise re- 
quire prohibitively long integrations). We find the follow- 
ing value for the differential Lindblad torque: ATlr m 
— 2.14Fo. We comment that, in order to evaluate Tq by 
means of Eq. ((62)) , we use a value of x s inferred from a 
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Table 2 

Parameters used in the different series of calculations. Each series consists of 21 runs, with n 
ranging from to 20. The viscosity or diffusivity therefore ranges from 2.10 - 9 a 2 f2 p to 
2.10 - 5 a 2 f2 p , in a geometric sequence. In the runs with an energy equation, whenever a 
parameters varies, the other is kept fixed with value 10 — 6 a 2 f2 p . This value was chosen to 
minimize the saturation due to either viscous or diffusive processes. 
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Figure 3. Total torque as a function of the time for the runs of set 
1 (isothermal equation of state). The dashed line shows an addi- 
tional inviscid run performed with the same set of parameters than 
the others. The viscosity ranges from 2 • 10~ 9 a 2 f2 p (lower oscilla- 
tory curves, darker) to 2- 10 — 5 a 2 O p (upper curves, lighter). On the 
electronic version the viscosity increases as the "wavelength" of the 
color decreases (the palette is approximately that of the rainbow). 

streamline analysis. We find xf° = 0.0153a, a value tha t 
is essentially independent of V (jMasset fc Ca soli 2009). 
The corresponding value of Tq is used to normalize the 
torque value in all our isothermal runs. 

5.3.2. Runs with an energy equation 

The results of set 2 and 3 are displayed respectively 
in Fig. [4] and [5j Calculations performed at low viscosity 
(therefore in set 2, see Fig. [4]) show a highly time varying 
behavior which prevents from extracting useful average 
values from these calculations. We shall discuss in a little 
more detail this behavior in section 16.5.11 It is straight- 
forward to estimate the differential Lindblad torque for 
these two sets: since there is no vortensity gradient, it 
suffices to perform a run with the same configuration, 
but with an isothermal equation of state. The gas being 
barotropic, the total torque amounts to the differential 
Lindblad torque, which we divide by 7 = 1.4 to obtain 
the differential Lindblad torque in the ru ns with an en- 
ergy equation (Ba ruteau &: Massetl [20081. We measure 
ATlr ss — 2.13Fo, where the value of x s that we adopt 
to estimate T is xf /^ 1 ' 4 . 

We make the following comment regarding these two 
sets: in an earlier version of these calculations, we noticed 
a long term drift of the torque value, which we attributed 
to a possible issue of boundary conditions associated with 
a relative drift of the disk's material with respect to the 
orbit. In order to get rid of this drift, we have adopted, 
for these two sets of calculations only, a profile of the 



Figure 4. Total torque as a function of time for the runs of the 
set 2 (varying viscosity). The color code is the same as in Fig. |3) a 
same color refers to the same value of n. For legibility purposes, the 
results of runs n = to n = 5 (darker curves, reddish colors on the 
electronic version) are averaged over a temporal window of width 
120 orbits, as they would otherwise display very large oscillations. 
For larger values of the viscosity a smooth behavior is observed. 
We note that the time averaging of the torque value distorts its 
initial behavior: without averaging, the torque of runs n = to 5 
would fall on top of the other curves for t < 50 orbits. 

kinematic viscosity that reads: 

v = vo (-) , (68) 

where Vq is the kinematic viscosity at the orbit. One 
can check that this prescription yields a vanishing radial 
drift of the disk's material for the surface density profile 
in r~ 3 / 2 adopted here. The runs performed with the 
prescription of Eq. (|68[) display a better behaved torque 
on the long term, which allows to measure an asymptotic 
torque value with a satisfactory accuracy. 

We note that the curves at low thermal diffusivity in 
figure [5] do not pile up at the value of the differential 
Lindblad torque, which is out of the range of the vertical 
axis. Although this does not contradict any first prin- 
ciple (since all these runs are performed with the large 
viscosity v = 10 _6 a 2 f2 p , so that they can sustain an un- 
saturated corotation torque), this is somehow surprising, 
as one would expect to be left with a vanishing corotation 
torque, when the thermal diffusivity tends to zero, since 
there is no initial vortensity gradient, and since the sin- 
gular vortensity created at the downstream separatrices 
tends to zero as the entropy is flattened over the horse- 
shoe region. This is not the case, however, and we will 
quantify this effect in section IB. 4. II 

5.4. Comparison with the reduced model 
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Figure 5. Total torque as a function of time for the runs of the 
set 3 (varying thermal diffusivity ) . The color code is the same as 
in Figs. [3] and [4] No temporal averaging is performed for these 
graphs. The dashed line shows the torque in the adiabatic case 
(vanishing thermal diffusivity). 

We have performed the same set of hydrodynamical 
calculations presented in section I5.2I in the framework 
of the reduced model detailed in section [31 We show 
hereafter how the results of this model compare to the 
results of the full hydrodynamical calculations. We have 
chosen x s to have the same value as the one we mea- 
sure for the runs of th e full hydrodynamical ca lculations, 
using the method of Casoli fe Massetl •: 200!);. which is 
x s sa 0.0153a. The mesh of the reduced model has 
N x = 4000 and N y = 100, and it extends from x = -0.3 
to x = 0.3, as in the hydrodynamical calculations. We 
add the estimate of the differential Lindblad torque given 
at section [5. 3. 1[ in order to compare directly the results 
of full hydrodynamical calculations and those of the re- 
duced model. 

We show in Fig.|S]the torque as a function of time in the 
barotropic case. Apart from the calculations which have 
a large viscosity, these results compare very well to the 
full hydrodynamical calculations (Fig. [3]), which shows 
that most of the long term torque behaviour observed 
in hydrodynamical calculations is accounted for by the 
horseshoe dynamics and a simple diffusion equation for 
the load. Even the slow decay of all torques with time, 
still noticeable at t = 500 orbits, is reproduced in the re- 
duced model. The agreement ceases when the viscosity 
is too large, namely for v > 2.10 6 . For these values, the 
timescale of the diffusion of load over the width of the 
horseshoe region, which is t v ~ x 2 s /v, amounts to ~ 20 
orbits or less, that is to say, in order of magnitude, the 
time required to perform a horseshoe U-turn or to estab- 
lish the horseshoe drag, as can be seen in Fig. [3] For these 
large values of the viscosity, the vortensity is therefore 
not conserved during a horseshoe U-turn, a property that 
is not captured by the reduced model, which assumes 
that the U-turns are instantaneously performed, so that 
they strictly conserve vortensity, by virtue of Eqs. (|38l) 
and ([39]). 

5.4.1. Cutoff function at large dissipation 

The discrepancy at large viscosity underlined above 
is already noticeable from the results over the first 
half libration time, i.e. over the first ~ 50 or- 
bits. One should indeed recover the linear estimate 
of the corotation torque in the lim it of large diffusion 
(|Paardekooper fc Papaloizoull2008f ). We first entertain 



Figure 6. Torque as a function of time for the 21 calculations of 
the run set 1, performed in the framework of the reduced model. 

the barotropic case, in the limit of large viscosity. 

In order to model in a simple manner this decay, we 
make the simple assumption that for a given viscosity, 
the unsaturated corotation torque is a blend of the un- 
saturated horseshoe drag as given by the reduced model, 
and of the linear corotation torque, so that: 

Tcr = £br H s + (1 - £fc)rcR,iin, (69) 

where £(, € [0, 1] is a blending coefficient. We can extract 
from our calculations the dependency of the £& coefficient 
upon the viscosity, provided we have an estimate of the 
linear corotation torque. The latter can be estimated 
by subtracting the differential Lindblad torque from the 
total torque at t rj 2 orbits, which yields TcRjin ~ l.OFn. 
We expect that for a sufficiently low viscosity, « 1, 
while Eb should tend to at large viscosity. This suggests 
the following form for £{,: 



1 + -ftTTu-turn/TV' 

where K = 0(1) is a dimensionless coefficient, ru-turn 
is the time scale for a horseshoe U-turn, and t v = x 2 s jv 
is the viscous timescale across the horseshoe region. We 
measure ru-turn from the run i = of the set 1, by eval- 
uating the time it takes for a fluid element to go from 
x = x s ,cf> = 1 rad to x = —x s ,4> = 1 rad. We find 
ru-turn ~ 1-1 • 10 2 ftp 1 - We compare the blending co- 
efficient given by the simulations to the trend given by 
Eq. ([70]) in Fig. [3 We see in figured] how the correction 
of the reduced model yields results which are in close 
agreement to those of the full hydrodynamical calcula- 
tions. From now on, we shall consider that the results 
of the full hydrodynamical simulations are correctly re- 
produced by the reduced model to which we apply the 
cutoff provided by the blending coefficient of Eq. (f70|) . 
This cutoff can be obtained from short runs, that are 
limited to the first half libration time after the planet is 
introduced in the disk. 

We repeat a similar analysis for the runs with an en- 
ergy equation. The decay of the torque at large dissi- 
pation can occur cither for a large viscosity or a large 
thermal diffusivity. Assuming the behavior to be separa- 
ble in viscosity and thermal diffusivity, we find from the 
analysis of the torque over the first 50 orbits that the fol- 
lowing ad hoc functions render account of the torque cut 
off, respectively at large viscosity and at large thermal 
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Figure 7. Blending coefficient of the isothermal series, as a func- 
tion of v isco sity. The solid curve is extracted from simulation s us - 
ing Eq. H69I I. whereas the dotted line shows the result of Eq. ( 1 701) ■ 
with K = 1. 




Figure 8. Tota l tor que inferred from the reduced model, to which 
the cutoff of Eq. d 70 1) is applied at large viscosity. This figure should 
be compared to Fig. [3] 



diffusivity: 



1 



1 + (itYu^turn/T,,) 2 ' 
1 

1 + Kru- tuin /T K ' 



with K = 0.2, (71) 
with K = 0.5. (72) 



6. SATURATION FUNCTIONS 

In this section we work out suitable analytic expres- 
sions for the perturbed vortensity and entropy fields at 
large time, and the corresponding asymptotic horseshoe 
drag estimates. We perform our analysis in the frame- 
work of the assumptions mentioned in section 13.21 to- 
gether with the boundary conditions specified in sec- 
tion 13.31 and we assume that the flow is steady state 
in the frame corotating with the planet. 

6.1. Timescale for the asymptotic regime 

In a low diffusion case, it can take a very long time for 
the vortensity or entropy to relax towards a steady state. 
For the lowest values considered here (v — 2A0~ 9 a 2 ilp 
or k = 2. 10~ 9 a 2 fl p ), the time required to relax towards 
a steady state over the whole mesh is: 



^rcla 



\R 



max, mm 



min(i^, k) 



10 7 orbits. 



(73) 



This time, which is longer than the lifetime of the disk if 
a > 1 AU, is much longer than the time scale accessible 
to fully hydrodynamical calculations, and is still beyond 
the reach of reduced models (we are able to run reduced 
models at most over ~ 10 5 orbits). In order to assess the 
validity of our analytic expressions, we therefore proceed 
in a two step manner: 

• We compare the results of the reduced model with 
the full hydrodynamical calculations at the final 
date reached by the latter, in order to assess the 
accuracy of the results of the reduced model. 

• We compare our analytic expressions with the re- 
sults of the reduced model at much larger times. 

6.2. Saturation in the barotropic case 



We note that Eq. (|34|) . which describes the evolution 
of the load in the barotropic case, is formally identical to 
Eq. (|36p , which describes the evolution of the entropy in 
the baroclinic case (one has to substitute the load with 
the entropy, the viscosity with the thermal diffusivity, 
and the gradient of load with the gradient of entropy). 
This analysis therefore also provides the information that 
we shall later need to evaluate the production of load at 
the downstream separatrix in the baroclinic case, as it 
requires the value of entropy at the separatrices. For the 
sake of definiteness, we consider the load in the remaining 
of this section. Adapting the results to the entropy is 
straightforward . 

6.2.1. Reduction to a convolution equation at low diffusion 

If we perform the change of variable u = x 3 ' 2 , Eq. (|34|) 
can be recast, for x > 0, as: 



8A _ d^L v_dL_ 
~9 v d^ + 3^~chl- 



(74) 



If we assume that L varies more rapidly than u, Eq. (|74l) 
can be recast as: 



8A 
~9~ 



d u L = v 



d 2 L 

du 2 



(75) 



Eq. ([75)1 is a diffusion equation in which y plays the role 
of the time (we note that A < 0, and that y decreases 
as time evolves for a fluid particle with x > 0). We 
denote with the value of the reduced load at y = 
27ra(l — £). We therefore have Lq = L R and L\ = L F . 
For a given initial configuration Lq in y = 2ira, the final 
configuration, after the fluid has executed a full orbit in 
the corotating frame, is given by: 



L x = L{y = 0) = L *G, 



(76) 



where * denotes the convolution product and where the 
Green's function G is given by: 



G(u) 



1 



ny/9va/\A\ 



exp[-u 2 |A|/(9i/7ra)] 



(77) 



The above formulation is valid when the width of the 
convolution kernel is small compared to the width of 
the horseshoe region (otherwise edge effects may be im- 
portant, and a significant diffusion from the region with 
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opposite of the field represented by the solid thick line 
(rear field) for x < x s , while they coincide for x > x s . 
We focus our analysis on the properties of the field in 
the vicinity of the separatrix. We aim at providing a 
universal function that represents the load in normalized 
coordinates in the low diffusion limit. In steady state, 
the field verifies the following equalities: 



L F (X)- 



■■ L R -k G 



L R (x) if x > x s ; 
—L R (x) otherwise. 



(83) 
(84) 



For convenience, we define the variable 



Figure 9. Perturbed load for different values of the azimuth. The 
horseshoe separatrix is at x s = 0.0153. The thick solid line, which 
has a discontinuity in x = x s (X = 0), represents the field just 
after the outward horseshoe U-turn (at the rear of the planet). 
The other curves represent the field value at different values of the 
azimuth (in geometric sequence). The thick dotted line shows the 
front load (at y = 0, just prior to the inward horseshoe U-turn.) 
It coincides with the rear distribution (solid line) for x > x a . The 
dash-dotted line shows the analytic expectation of Eq. J 103 tt ■ The 
smooth solid line represents the azimuthal average of the load. 

x < 0, which drifts towards increasing values of y, may 
occur). This condition translates into: 



X = 



9van 

1M 



< x° a 



We define: 



av 



(78) 



(79) 



which represents, within a numerical factor, the ratio of 
the libration and viscous diffusion time scales across the 
horseshoe region. The condition of Eq. (1781 therefore 
translates into: 

z v < 0.026 = z{ (80) 

Similarly we define 



OK 



(81) 



which characterizes the diffusion of entropy. We say that 
a given value of viscosity (or thermal diffusivity) falls in 
the regime of low diffusion when Eq. (|80|) is satisfied (or 
when z K < 0.026). This means that a localized feature is 
spread over less than the horseshoe zone width after half 
a libration time. When this condition is satisfied, we can 
disregard the shear and write the final distribution of the 
load, near the separatrix and after a synodic period (half 
a libration), as a convolution product directly in x, so 
that the Green's kernel reads: 



G(x) 



1 



V47W70 



exp(-z 2 /4i/To), 



(82) 



where To = 4:Tra/(3flpX s ) is half the horseshoe libration 
time near the separatrix. 

6.2.2. Properties of the stationary flow at low diffusion 

Figure |9] shows a set of radial cuts of the load for dif- 
ferent values of the azimuth, in a regime of low diffusion. 
One easily checks that this field obeys the boundary con- 
ditions given by Eqs. (J37J, ([3"9|). and (j44j), with A = 0: 
the field represented by the dotted line (front field) is the 



(85) 



Eqs. (|8"3")l and (1541 read, expressed in term of the variable 
X: 

L F = L R *K where K{X) = cxp(-X 2 )/V^ (86) 
and 

(87) 



L F {X) = sgn(X)L R (X). 

We consider the function /i that verifies Eqs. and 
(|87p and that obeys the following constraints: 

/i(0) = 1 
lim f\(X) = 0. 

X—^ — oo 



(89) 



This last condition comes from the fact that in the low 
diffusion limit, the load profile is essentially flattened 
over the horseshoe region, and therefore tends to zero 
towards corotation (i.e. for X —¥ — oo). This function 
can be obtained by iteration of the transforms T% and Ti 
respectively defined by: 



and 



Ti(f) : X ^ 



T 2 (f) 



f(X) 
~f(X) 

f*K 

(f*Ky(oy 



if X > 
otherwise 



(90) 



(91) 



where the division of the right hand side is a normaliza- 
tion imposed to ensure that Eq. (|88"1) be satisfied. The 
action of Ti represents the U-turns, whereas T 2 describes 
the viscous diffusion between two subsequent U-turns. 
We iterate these transforms starting from an arbitrary 
initial function /, so that iteration number n is denoted 

9n- 

9 n = T 2 oT l [T 2 oT 1 {f)]. (92) 



The iterates are found to converge so that we can adopt 



h 



lim r , 



fjn 



In this iterating procedure, we set 



g n to below an arbitrary large, negative value of X, 
in order to satisfy Eq. (|59"|). It is straightforward to 
check that the couple of functions L F (X) = fi(X) 
and L R (X) = sgn(X)fi(X) thus constructed satisfies 
Eqs. (US])- dHHJ) - The function fx(X) is depicted in fig- 
ure 1101 The actual load distribution, and ultimately the 
horseshoe drag, can be obtained by a proper scaling of 
the function f\. In order to perform this scaling, we need 
an additional relationship on the radial derivative of the 
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Figure 10. Normalized load distribution in the low diffusion limit. 
The solid line shows /i, while the dashed line shows its derivative. 
The horizontal dotted lines show the value of f\ (0) and the asymp- 
totic value of f'(X) for X —> oo. 

load. Integrating Eq. in y over the range [0,2™], 
and using the boundary condition given in Eq. (1371) . we 
note that: 

p2na 

/ L"dy = for any |sc| > x s , (93) 
Jo 

that is to say the azimuthally averaged slope of the field 
is uniform outside of the horseshoe region (even though 
the timescale required to reach this state may be long, as 
underlined in section RTTI ) Owing to the boundary condi- 
tions given at Eqs. (f4"2"j) - p5)l . this azimuthally averaged 
slope is therefore necessarily: 

<£') = (94) 

where the symbol (. . .) denotes the azimuthal average. 
We can estimate the azimuthally averaged slope of our 
normalized function in X = (i.e. at the separatrix) in 
two different, equivalent manners: 

• The function fi(X) represents the load distribu- 
tion in front of the planet (y = 0). It does 
not intrinsically contain sufficient information to 
yield the azimuthally averaged slope. If we note 
£ = 1 — yj (2ira) a normalized value of the azimuth 
(which is at the rear of the planet and 1 in front 
of the planet), then, in normalized coordinates, the 
load distribution f^(X) at azimuth £ is given by: 

f( =T 1 (f 1 )*K ( = f *K i , (95) 

where 

iQ = ^L=exp(-X 2 /0- (96) 
The azimuthally averaged slope is then given by: 

(4(°))«= / pi(/i)*jy'(o)de (97) 

Jo 

We note that when £ — > 0, the slope can be arbi- 
trarily large, and that the function / = is 
discontinuous in X = 0, so that a direct estimate of 
the integral of Eq. ([97]) may be inaccurate. Rather, 
if we note jq — /i(0) « 0.52, then we decompose 
/o as follows: 

MX) = sgn(X) 7o + [h(X) - 7o]sgn(X). (98) 



The first term of the right hand side of Eq. is 
a step function, whose contribution to Eq. (|9T[) can 
be evaluated analytically, while the second term 
of Eq. (j98| is continuous in X — 0, so that an 
accurate estimate of its contribution can be ob- 
tained by direct summation of tabulated estimates 
of [(fi(X) - 7o )sgn(X)] * K^. This yields the fol- 
lowing estimate of the azimuthally averaged slope, 
that we denote with 7 i: 

7i = (/e(0)>? « !-48 (99) 

• One can notice that the slope of the load at any az- 
imuth tends towards a unique value at large X (and 
therefore towards the azimuthal average). This can 
be inferred from examination of Fig. [§J where one 
can see that the different curves cluster at large, 
positive values of X. The asymptotic value of 
f[{X) at large X is therefore also 71. This is rep- 
resented in Fig. QUI in which we see directly that 
71 « 1.48. 

The load at the upstream separatrix is scaled as follows. 
It reads: 

^w-^fe)- (100 » 

where fj, is the constant that we aim at estimating. The 
radial derivative of the load at the separatrix reads, using 
Eqs. ([Mf and ([9"9|: 

d x L F \ x = —L' 00 (101) 
on the one hand, and, using Eqs. {55} and ([55]): 

d x L F \ =^= (102) 

on the other hand. Using Eqs. (|101[) and (|102[) we infer: 
L F {x) = ^V^oh (^^) ■ (103) 

Figure [9] shows this estimates with a dot-dashed line. It 
compares reasonably well to the true load distribution, 
which has z v ~ 0.014 in the example shown (half the 
threshold value of the low diffusion regime). Of partic- 
ular interest is the value L s at the separatrix, upstream 
of the U-turn. Wc have: 

L s = L F {x s ) = M /i(0) = L'J^^UFq. (104) 

7i 

Using Eq. (|43[) , Eq. (|104|) can be transformed using a for- 
mal analogy to give the entropy upstream of the horse- 
shoe U-turns, at the separatrices: 

S s = S F {x s ) = S'^V^Fo. (105) 
7i 

This is a key quantity to evaluate the singular production 
of load at the downstream separatrices, given by Eq. (l48|) . 
In the low thermal diffusion regime, the entropy related 
torque therefore scales with k 1 / 2 . 
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Figure 11. Load value on t he separatrix, at y = 0. The dashed 
line corresponds to Eq. (1 10411 up to values of viscosity that yield 
the value of the load in the unperturbed disk (i.e. for z„ < z%), 
and to the value of the unperturbed disk for z v > zS. The solid 
line corresponds to the results of the reduced model. We see that 
the analytic expression accurately accounts for the results of the 
reduced model in the low diffusion regime [z v < zj). For larger 
viscosity, some discrepancy is observed, with a maximal relative 
error of about 20 %. 

6.2.3. Case of the large diffusion regime 

In the low diffusion regime, the value of the load at the 
upstream separatrix given by Eq. (|104|) is much smaller 
than its value in the unperturbed disk, which is 



L_ — L„ 



(106) 



This is due to the fact that the advection of load within 
the horseshoe region tends to flatten its profile over the 
horseshoe region. Reciprocally we expect Eq. (I104[) to 
break down when it yields values larger than the unper- 
turbed value, which happens for 



> 



37i 2 

167T 7 2 



0.48 = z% 



(107) 



Above this value we expect the diffusion to be sufficiently 
strong to impose at any instant in time the unperturbed 
load profile, so we expect the load at the separatrix to 
saturate at the value LP S . Figure [TT1 shows how this ex- 
pectation compares with the results obtained from the 
reduced model. We see in this figure that the dashed 
line provides a reasonable approximation to the actual 
value of the reduced model. We shall adopt henceforth 
the following approximation of the load value at the sep- 
aratrix, upstream of the U-turn: 



L s = L'oo min ( x s , —\fAi>T 
V 7i 

for the load in the barotropic case, and 



S s = S'ec min ( x s , —s/Ikto 

7i 



(108) 



(109) 



for the entropy value at the separatrix upstream of the 
U-turn, in the case with an energy equation. 

6.3. Saturated horseshoe drag 

Integrating Eq. (IMl) from y — to y — 2ira we obtain, 
in a steady state flow: 

2Ax[L R (x) - L F (x)} = 2naisdl 2 (L), (110) 



which yields for any |x| < x s , using Eq. (|45l) : 

- AAxL F {x) = 2navdl 2 {L). (Ill) 

Integrating by parts, and using the symmetry property 
of Eq. ((44]) . this allows to transform the horseshoe drag 
expression of Eq. ([59]) into: 

V = %TrB^ V a 2 [x s d x {L) s -{L s )], (112) 

where we have made use of the fact that (L) x= q — 
0, which arises from the sym metry property given by 
Eq. (|gjl. As in iMassetl (|2001[ ). we find that the horse- 
shoe drag reduces to a term to be evaluated at the sep- 
aratrix. The first term of the bracket of Eq. (|112[) can 
be expressed using Eq. (|94|) . The second term requires 
the estimate of the azimuthal average of the load at the 
separatrix, much in the same way as we evaluated the 
azimuthal average of its derivative in section 16.2.21 We 
obtain: 

(A(0)) e « 0.36 = 7270 with 72 » 0.68, (113) 
hence we have, in the low viscosity limit, using Eq. (|104p : 



727o 
7i 



(114) 



Using Eqs. ©, (JEJJ), (J94J) , ([TT2]) and (fTT4]) we can write 
the horseshoe drag, in the low diffusion limit, as: 



r = 8nBvaVTlQX s 



47270 



: 2 / nQ, p vaVY,QX s {\ 



7i 



1/2 



.1/2 



(115) 
(116) 

where we have restricted ourselves to the Keplerian case, 
and where we have used the property: 



47270 
7i 



1/2 



0.99 w 1. 



(117) 



At very low viscosity, z u — ¥ and r oc v. In the low 
diffusion regime, Eq. (|116[) can be recast as: 



8tt 



vroz^i-zy 2 ), 



(118) 



where To is given by Eq. (|62j) . Figure [12] shows how 
Eq. (|118[) compares to the actual horseshoe drag given by 
the calculations in the framework of the reduced model. 
This expression is bound to fail at some point since it 
yields a negative result for z v > 1. It has a maximum 
at z v = 4/9, which is (327r/81)Vr ~ 1.24r . We see 
that the agreement between the calculations and this ex- 
pression is remarkably good up to this value of z (which 
roughly coincides with z%)- Figure [T^] shows that be- 
yond z%, the flat approximation provides a good approx- 
imation to the horseshoe drag, with a relative error of 
about ~ 15 % at most, for realistic values of the viscos- 
ity. We therefore adopt the following expression for the 
barotropic horseshoe drag: 



87T 



VY z v F(z v ), 



where 



F(X) 



{ 1-X 1 / 2 
\ 4/(27X) 



if X < 4/9 
otherwise. 



(119) 



(120) 
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Figure 12. Horseshoe drag as a function of viscosity, normalized 
to VTo. The solid line shows the results of the reduced models. 
Note that the viscosity range has been significantly extended to- 
wards larger viscosities in order to capture the decay of the horse- 
shoe drag past the peak around z u ~ 3 and its convergence to- 
wards t he s tandard unsaturated barotropic horseshoe drag given 
by Eq. J 6 IP - The upper dotted curve shows the linear expression 
(87r/3)Vroz, while the lower dotted curve — w hich coincides with 
the thick dashed curve for z < z% — shows Eq. 11181 . The two ver- 
tical dot-dashed lines show the critical values of the viscosity that 
correspond to zf and z^. The thick dashed curve shows Eq. ||119| I 

We note that, in the case of any field that obeys a dynam- 
ics similar to that of the load in the barotropic case, i.e. 
that obeys the governing equation and boundary condi- 
tions given at sections 13.21 and 13.31 (without the produc- 
tion of singular stripes at the downstream separatrices), 
we have the following relationship: 

x s (d x L) s - (L) s = x.L'^Ffa) (121) 

6.3.1. Comparison with hydrodynamical calculations and 
previous results 

We conclude this section with a comparison of the re- 
sults of hydrodynamical calculations with the analytic 
expres sion, and with the expression given by iMassetl 
(|200lD . Figure [13] displays the results of full hydro- 
dynamical calculations at a relatively early stage (t = 
500 orbits), as well as the results of the reduced model, 
at the same date and at a much larger time. Some differ- 
ence can be noted between these two dates, which shows 
that the corotation torque at t = 500 has not reached a 
steady value, in agreement with the slight decay notice- 
able in the torque curves of figure [3] at large time. The 
agreement between the full hydrodynamical calculation 
and the reduced model at the same date is excellent. We 
note in passing that a slight (10 %) rescaling has been 
applied to the viscosity value of the reduced model. The 
saturation state of the horseshoe region depends indeed 
upon the ratio of the libration timescale to the viscous 
timescale. This ratio is slightly larger in a real situation 
than in the corresponding case for the reduced model. 
Indeed the horseshoe streamlines are not located at a 
constant distance from corotation. The existence of the 
potential's indirect term renders the streamline slightly 
more narrow near opposition (0 = n). This effect is more 
important in thicker disks, and is marginally relevant for 
the disk aspect ratio of h = 0.05 adopted in our calcula- 
tions. This has two consequences: 

• the libration timescale is slightly larger, as the ma- 
terial drifts slower near the opposition. 
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Figure 13. Total specific torque as a function of viscosity. The 
upper solid line with diamonds shows the torque at t = 500 or- 
bits in full hydrodynamical calculations. The dashed and dash- 
dotted lines show the results of the reduced model, respectively 
at t = 500 orbits and t ~ 10 5 orbits. The lower solid line shows 
the analytic expectation given by Eq. 111911. w hile the dotted line 
shows the analytic pre dict ion of Massct ( 2001]). The torque cutoff 
function given by Eq. 1701 1 is applied to the reduced model and to 
the analytic expression, except that of the dotted line. 

• the viscous timescale is on the average slightly 
shorter, as the azimuthally averaged horseshoe 
width is smaller. 

Both effects act at increasing the role of diffusion. A 
streamline analysis of a typical output of the hydrody- 
namical calculation shows that the libration to viscous 
diffusion timescale is ~ 1.1 larger than the crude esti- 
mate of the reduced model. Hence the plots of the re- 
duced model results and of the analytical expectations of 
figure f!3l are performed using a value of the viscosity 1.1 
times larger than the current viscosity, which amounts 
to slightly shifting leftward the corresponding graphs on 
the figure. We make the following additional comments: 

• The difference between the analytical expression 
and the hydrodynamical calculations is entirely ac- 
counted for by the date effect: reaching a full 
steady state is beyond the possibilities of hydro- 
dynamical calculations, and the asymptotic torque 
is slightly more saturated than at t = 500 orbits. 

• The expression of Masset (200l|), even if it captures 
correctly the range of viscosities over which a tran- 
sition is observed from full saturation to desatura- 
tion, does not quite provid e as a ccurate an estimate 
of the torque value as Eq. Q119p . This estimate was 
based upon the simplifying assumption that the 
material distribution within the horseshoe region 
has an axisymmetric distribution, which leads to an 
underestimate of the flow of load on the horseshoe 
U-turns. In addition, Eq. (11191) is much simpler to 
evaluate than is Eq. (9) of lMassetl (|2001l ) , which in- 
volves the evaluation of Airy functions. Since one 
of the main purposes of the present work is to pro- 
vide torque expressions that can be used in studies 
of planetary population synthesis, it is desirable to 
use an expression which minimizes the computa- 
tional cost of the torque estimate. 

We also notice that the steady state horseshoe drag can 
be larger than the unsaturated horseshoe drag, given by 
Eq. (|6ip. This is apparent both for the reduced model 
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Figure 14. Load profile in front of the planet (left) and corre- 
sponding map of the load (right). In the left plot, the grey shaded 
area shows the actual load profile that imposes the horseshoe drag, 
while the hashed one shows the load profile that corres pond s to 
the standard, unsaturated horseshoe drag given by Eq. 16H . In 
the right plot, the dark arrows sketch the motion of fluid, and 
the light arrows sketch how viscous diffusion spreads the load out 
of the horseshoe region (delimited by the dashed line) and to the 
upstream region. 

and for the full hydrodynamical calculations. This be- 
havior can be qualitatively understood as follows: when 
the viscosity is large, so that the diffusion timescale 
across the horseshoe region is shorter than the libration 
timescale, the flow downstream of the U-turns can act 
back on the upstream flow by diffusion, as sketched on 
the right plot of figure In a steady state situation, 
this can result in a load profile in front of the planet 
significantly above the unperturbed profile, as depicted 
in the left graph of figure [TU This behavior also yields 
a value of the load at the separatrix slightly above the 
unperturbed one, as can be noticed in figure QTJ One 
can estimate the excess of horseshoe drag in the situa- 
tion depicted in figure [TJ] by assuming that the load over 
the upstream part of the horseshoe region (correspond- 
ing to the grey shaded area) is uniform and equal to the 
value that it has at the separatrix in the unperturbed 
disk. Using Eq. ([59|) . it is straightforward to see that the 
horseshoe drag then supersedes the standard value by a 
factor 4/3. This is in good agreement with the peak value 
of the horseshoe drag in figure [T2] We comment that this 
excess is less pronounced in full hydrodynamical calcu- 
lations, since it corresponds to values of the viscosity for 
which some cutoff already occurs. 

6.4. Saturation of the entropy related torque 

We examine in this section the saturation properties of 
the entropy related torque, and exhibit an analytic ex- 
pression for the steady state entropy related corotation 
torque. Describing the saturation properties of the en- 
tropy related torque amounts to finding the steady state 
pattern of load that settles over a long time scale, when 
it receives a singular kick described by Eqs. (|3"5|) and (|39[) 
at each horseshoe U-turn, much as in the barotropic case. 
As explained in section 15721 we split the general case into 
a bulk term that has a non-vanishing slope at infinity, 
and an edge term that has a flat profile at infinity but 
undergoes the creation of singularities at the downstream 
separatrices. We focus firstly on the bulk term, and we 
disregard the singularities. 

6.4.1. Vortensity-Entropy viscous coupling 

For a barotropic disk we have made the simplifying as- 
sumption that the load (or vortensity) evolution could be 
described by a diffusion equation, i.e. we have neglected 
the last term of Eq. (|2~51) , so that the evolution equation 
of the load was described by Eq. (|2"6")) . Nevertheless, the 
last term of Eq. (|25[) may be arbitrarily large at a con- 
tact discontinuity, and we may expect this term to play 



a significant role in runs with an energy equation, at low 
diffusivity, when the contact discontinuity does not sig- 
nificantly decay along the separatrices. Under these cir- 
cumstances, the contact discontinuity may act as a source 
of load, which ultimately exerts a corotation torque. In 
order to work out the contribution to the torque of the 
load thus created, we make two amendments to the calcu- 
lations developed in the previous section: (i) the coupled 
evolution of the load and of the entropy is different, and 
described by the evolution of an eigenfunction that is a 
linear combination of load and entropy; (ii) the integra- 
tion by parts presented in section 16.31 must be modified 
accordingly, and yields a torque expression different from 
that of Eq. (|112|) . We briefly describe hereafter the cor- 
responding calculations. 

The evolution of the system, which is governed by 
Eq. (I3"5j) and Eq. (|3"6"1) . can also be described using the 
eigenfunction U defined by: 



U = L 



70 - «) 



(122) 



By virtue of the radial boundary conditions of Eqs. (|42|) 
and (|43|) . the asymptotic slope of U at large |x| is given 
by: 

lim ^Z7= Z7^ = ^ . (123) 

x^±oo a v — k a 

The governing equations of U are therefore similar to 
those of L in the barotropic case, except that U has a 
different slope at infinity. Namely, the evolution of the 
system is described by: 



d t U 



2Axd y U 



T Til 

vU 



(124) 



together with Eq. (|36|) and the boundary conditions of 
Eqs. (|123p and (j43|) . Also the function U undergoes the 
creation of singularities ±Lo at the downstream separa- 
trices, as described in Eqs. fl38p and (f3l?|). since L has 
a unitary weight in Eq. (|122p . As explained above, we 
disregard these singularities in this section, and we fo- 
cus on the bulk term exclusively. The contribution of 
the singular creation of load will be analyzed in the next 
section. 

We integrate Eq. (|55"j) from y = to y = 2na, assuming 
a steady state. This yields: 

- AAxL F = 2-KaviL") - — <S"). (125) 

7 

We multiply Eq. (|125l) by 8B 2 ax , and integrate it by 
parts between and x s : 



HS 



x s (U' s ) ~{U.) + 



(x s (S' s ) - (S s )) 
(126) 



Using Eq. (|121l) . which applies to a field such as U, as 
we outlined in section T6.31 we write: 



and 



x s (U' s )-{U s ) = U' 00 x s F{z v ) 



x s (S' s ) - (S s ) = S' 00 x s F(z K ). 



(127) 



(128) 
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Using Eqs. (jMJ), ([125]) , (TT2B)) . (TT2T)) and (p~^HI) we write: 



8tt 



HS 



VF{z v ) - 2S 



z v F{z v ) - z K F(z K ) 



(129) 

The last term in the bracket of Eq. (|129|) represents the 
additional contribution due to the vortensity-entropy vis- 
cous coupling in Eq. (|36|) . By continuity, in the particular 
case v = k (unitary Prandtl number), the torque is: 

r HS = y TqZ " \ vf{z ^ ~ 2Sd ^ ( 13 °) 



Eq. (|130[) constitutes the bulk horseshoe drag in the 
general case. This bulk term renders account accurately 
of the corotation torque that we obtain at low ther- 
mal diffusion and finite viscosity (runs set 3, see sec- 
tion 15.3. 2|) . In this situation, the entropy gradient virtu- 
ally vanishes across the horseshoe region, so the corota- 
tion torque essentially reduces to the bulk term, which 
reads, in this particular situation with z K j z v — ¥ and 
V = 0: 

Ths = yr (-25)F(z I/ ). (131) 

The torque is therefore — 25/ V times that of the 
barotropic situation with v = lO~ 6 a 2 Q p , which is what 
one can infer from the Figs. [3] to [5] 

We now turn to the edge term which stems from the 
vortensity sheets at the downstream separatrices. We 
comment that in order to isolate its effect in numerical 
simulations, it was necessary to firstly establish Eq. ([1301) 
so as to subtract it from the horseshoe drag, since, in gen- 
eral, Eq. (|130j) shows that there is always a bulk term 
associated to the edge term (i.e. when there is a non- 
vanishing entropy gradient), even in disks with a flat 
vortensity profile. 

6.5. Corotation torque associated to the singular 
production of load 

As we have already considered in the previous section 
the bulk term in the general case, which stems from a 
non- vanishing large scale slope of U, we now consider 
the complementary situation and isolate the effects of the 
singular production of load at the downstream separatri- 
ces, and assume a vanishing large scale slope: = 0. 
Owing to the linearity of the flow governing equations, 
the total torque is the sum of both contributions. 

6.5.1. Considerations about the inviscid case 

The special case of an inviscid disk deserves some dis- 
cussion. It has already been noted in numerical simula- 
tions that in inviscid disks with thermal diffusion, the 
corotation torque tends to zero after several libration 
timescales, so that some v i scosity is needed to prevent 
saturation (iBaruteaul 120081 [Paardckoopcr & Papaloizou 
120081: Kiev fc Cridall2008[ ). This is expected on general 
grounds, as a sustained exchange of angular momentum 
between the horseshoe region and the rest of the disk, de- 
scribed by the Navier Stokes equations, necessarily relies 
on a non- vanis hing stress tensor, regardless of the th er- 
mal diffusivity ([Masset fc Casolill2009t [Massetil2008allbh . 
This simple remark provides a strong clue about the ori- 
gin of the entropy related torque. It had been initially 
argued that this component of the horseshoe drag came 




Figure 15. Perturbed entropy in a nearly inviscid disk with ther- 
mal diffusion. The grey level is the same for both plots. The upper 
plot corresponds to a situation during the first half libration time, 
i.e. to the fully unsaturated horseshoe drag, while the lower plot 
corresponds to approximately 12 full libration times. 

from the under- and over-dense regions that appear at 
the ends of the horseshoe region as the result of material 
conserving its entrop y whereas the pressure remains es- 
sentially unchanged flP aardckoo per fc Papaloizou! 120081 : 
iBaruteau fc Massetl 120081 ) . More recently, it has been 
argued that this effect does not contribute to the en- 
tropy related torque, which exclusively arises from the 
singul ar production of load at the downstream separa- 
trices (Massct & Casoli 2009). Figure [15] illustrates why 
this has to be the case indeed. It shows maps of per- 
turbed entropy corresponding to the run of set 2 with 
viscosity v — 2 • 10~ 9 a 2 Q p (so that it can be consid- 
ered, for our purpose, as an inviscid run). The thermal 
diffusivity in all the runs of this set is such that the ther- 
mal diffusion time is longer than the horseshoe U-turn 
time (so that material performing the U-turns behaves 
adiabatically) but slightly shorter than half the libration 
timescale, so that a fixed pattern of perturbed entropy 
subsists forever in the horseshoe region (associated with a 
corresponding pattern of perturbed density). Neverthe- 
less, the corotation torque tends to zero (on the average, 
since the situation is very messy owing to the presence 
of vortices), as it has been seen in figure [4] This rules 
out the possibility that this localized pattern contributes, 
even partially, to t he entropy related torque, as was re- 
cently suggested bv lPaardekooper et al.l (|2009f) . If it did, 
the corotation torque could not possibly saturate at null 
viscosity and finite thermal diffusion, which is in con- 
tradiction with first principles. With this in mind, it 
is instructive to undertake a simple minded analysis of 
the inviscid case. We consider a planet "switched on" 
in a initially unperturbed disk without vortensity gradi- 
ent. A singular amount of load ±Ao begins to drift along 
the downstream separatrices. After a time To, the singu- 
lar load has executed a complete orbit in the corotating 
frame. We can consider that the singular amount of load 
has spread radially as a vanishingly narrow Gaussian un- 
der the action of a vanishingly small viscosity. In an 
idealized situation, the half part of the Gaussian packet 
that is outside of the horseshoe region circulates in the 
outer (inner) disk, whereas the other half is shunt on a 
horseshoe U-turn. In practice, however, the horseshoe 
region is not strictly closed, especially at early times, 
and can have a sligh tly different width in front of the 
planet and behind it (|Casoli fc M assct 20091) . We there- 



20 



F. S. Masset & J. Casoli 



fore entertain a slightly more general situation in which 
a fraction 1/2 — x of the singular load circulates, whereas 
the remaining fraction 1/2 + x is shunt toward horseshoe 
U-turns, as depicted in figure OH We denote A„ the 
load at the front separatrix, downstream of the U-turn, 
over the time interval [utq, (n + l)ro[. The load value at 
this location changes every To, and so does the torque 
by virtue of Eq. (|39|) . Figure [TBI shows the relationship 
between A„_i and A„: 

A„ = (-- X j A„_i-f - +x) A„_i+A = -2xA n _i+A . 



From this recurrence relation we infer 

ifX 



[1 - (-2 X )" +1 ] if x ^ -1/2 
-1/2, 



A„ = t i+2x 

(n + 1)A 



(132) 
(133) 



and the corresponding torque value is, using Eqs. (|59p 
and fT33]l : 



r„ = W\A\B 2 ax 2 s Ao(-2 x y 



(134) 



In the general case (|x| < 1/2), the torque tends to zero 
at large time. It does so after just half a libration time 
in the particular case x = 0. If x > 0, it oscillates about 
0, whereas it tends monotonically to if x < 0- The 
ideal situation (x = 0, i.e. the load is equally split at 
the stagnation point between circulation and horseshoe 
U-turn) is more likely to be found in a "quiet" situa- 
tion, for instance at large softening length. We expect 
in that case that the torque quickly saturates (note that 
in practice, the value of Ao changes over time until the 
entropy field reaches a steady state, so that the thermal 
diffusion timescale determines the minimum amount of 
time needed to reach the torque saturation). 
The peculiar cases |x| = 1/2 are also of interest. 

• If X = 1/2 (which means that the whole amount of 
load executes horseshoe U-turns), the torque indef- 
initely oscillates about 0. The field does not con- 
verge towards a steady state. However, the time 
averaged torque does converge to zero, which can 
be considered as a weak version of the torque sat- 
uration. 

• If x = — 1/2 (which means that the load produced 
at the stagnation point never undergoes horseshoe 
U-turns), the load increases (decreases) steadily at 
the inner (outer separatrix). This rather artifi- 
cial case therefore corresponds to a steady flow of 
material from one side of the horseshoe region to 
the other side, which gives rise to a non vanishing 
corotation torque, much as in type III migration 
(|Masset fc Papaloizoull2003[ ). 

Therefore, apart from the artificial case x = — 1/2, the 
torque is bound to saturate. The manner in which it does 
so depends on the value of Xi hence it depends on very 
local properties of the flow in the vicinity of the stag- 
nation point. It can either oscillate significantly before 
saturating, or it can tend to zero quickly and monoton- 
ically. This simple minded analysis therefore illustrates 
the fact that, even if the asymptotic value does not de- 
pend on the details of the flow, the manner in which the 
torque saturates is highly sensitive to the flow topology, 
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Figure 16. Sketch of the load at a given instant in time in an 
inviscid disk. The x-axis is the azimuth and the y-axis is the radius. 
The horizontal solid lines show the separatrices, while the dashed 
line shows the corotation. 

and it may differ from one particular case to the other. 
This is quite in contrast with the behavior of the bulk 
term, in which the torque value is dominated by phase 
mixing, and has a temporal behavior much more robust 
and predictable. 

6.5.2. Universal form of the load in the low viscosity limit 

We seek the steady state distribution of the load in the 
low diffusion limit. This amounts to finding a distribu- 
tion which has a vanishing slope for \X\ oo and which 
satisfies Eq. (|83[) and the following relationship: 



L R (x) = -L 5(x-x s ) + 



L F (x) if \x\ > x s : 
—L F (x) otherwise. 



(135) 



Using the variable change of Eq. 
conditions as Eq. ([86]) and: 



we can recast these 



L R (X)= S gn(X)L F (X)+Z6(X), 



(136) 



where £ = —Lq/A. We seek a couple of functions 
(L F , L R ) that verifies Eqs. (JHSD and ([Hp) with £ = 1. It 
is straightforward to check that if one adopts for L F {X) 
the function 9l (X) defined by 9l {X) = f{(X)/(2>y ), 
where /i is the function that has been obtained in sec- 
tion [6X2l and for L R {X) the function given by Eq. (fl36l) . 
then the functions L F and L R so constructed verify 
Eq. (H|. The function f[(X) is represented in Fig. [10] 
For an arbitrary value of Lq, we have therefore: 



L F (X)=^ gi (X) = -^ 9l (X). 
The horseshoe drag of Eq. ([59]) can be written 



Ths = 2 / x L(x)dx + L x s , 
'a 



(137) 



(138) 



which we expand, to first order in y/4i>To/x s , as 

Ths 
8\A\BY, a 



/0 rO 
L F (X)dX-AAL x s / X gi (X)dX, 
- oo J — oo 

(139) 



A property of gi(X) is that: 

9l (X)dX = / 1 (0)/(2 7 o) = l/2, 



(140) 
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so that the horseshoe drag expression above reduces to: 

,o 

r HS = -32\A\BZ aV^L x s / X 9l (X)dX. (141) 

J — oo 



The integral on the right hand side of Eq. (|141[) is ap- 
proximately equal to —0.176, hence we can reduce the 
torque expression to, using Eq. (|47|) : 



2.82 



(142) 



The above expression, when compared to the results 
of the reduced model (not reproduced here), correctly 
accounts for a v 1 / 2 dependence, but overestimates the 
horseshoe drag of the reduced model by a significant fac- 
tor (about 1.6). The reason for this can be tracked back 
to the simplification of Eq. f| 139[) arising from Eq. (|140[) , 
and it illustrates again the high sensitivity of the torque 
edge term to the detail of the processes that affect the 
load advected along the separatrix. Whereas Eq. (|140|) 
is exact in the framework of the set of Eqs. (|136[) and 
, which come themselves from Eq. (jT5j) , it is only ap- 
proximate when the shear is retained, which means that 
a singular amount of load drifting along the separatrix 
is not exactly split into two equal fractions on the inside 
and on the outside of the separatrix. It is straightfor- 
ward to estimate the magnitude of this effect. We define 
the true load L t (X) distribution as the approximate one 
determined above plus a residue SL(X) that we want to 
characterize: 




L t (X) = L(X) + 6L(X), 



(143) 



and we make the assumption that |<5L| <C \L\. If we inject 
this into Eq. (|74|) , we obtain, after integrating from y = 
to y = 2na: 



Figure 17. Entropy related torque (edge term), as a function of 
viscosity. The solid line shows the value inferred from the runs 
of se t 2 (to which the estimated bulk horseshoe drag given by 
Eq. (1 1 29 1) has been s ubtra cted), while the dashed line shows the 
results of Eqs. I ll50t - ||151l l. The value chosen for TcRlm is 1*1 /3, 
given by estimates of the torque value only 2 orbits after the inser- 
tion of the planet in the disk. The inset plot represents the same 
quantity with a linear scale on the j/-axis. 

which agrees with the results of the reduced model within 
20 %. 

6.5.3. Extension to arbitrary viscosity 

Eq. (|149p remains valid as long as the viscosity is suf- 
ficiently small. Results obtained in the framework of the 
reduced model (not reproduced here) show essentially 
two branches: the branch in 2/ 1 ' 2 , corresponding to the 
regime of low viscosity, and a flat branch at large vis- 
cosity. These two branches intersect at z v — 0.3, hence 
we adopt the following expression for the horseshoe drag 
edge term: 



r — r 2 £„ + TcR,lin\\ — 



(150) 



^\A\6L F = 2ira-d u (L). 

6 U 

The horseshoe drag of Eq. (|139|) then becomes: 



(144) 



where 



T 2 = 2.16Vmin(,z 1/ ,0.3)ri 



r 



HS 



i\A\BE a 



/o i-O 
SL F (X)dX-4AL x s / X 9l (X)dX. 
-oo J — oo 



(145) 

In the low diffusion regime, the first integral of the right 
hand side of Eq. (|145|) reduces to: 



f 



SL F (X)dX 



A\A\x] 



(146) 



Using the relationship (L) s — 71 L F , we end up with an 
estimate of the additional contribution to the horseshoe 
drag due to 8L (the first term of the right hand side of 
Eq. (11451) : 

§Tji S = 16\ A\B-E ax 2 s f ' SL F (X)dX (147) 
Jo 
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Adding up Eqs. (|142|) and (|148p . we eventually get 



(148) 



HS 



2.ny|zy 2 r 1 «2.i64/ 2 r 1 , 



(149) 



(151) 

represents the analytic trend derived in the previous sec- 
tion, saturated at z v — 0.3, and where e v , given by 
Eq. (fTTl) , is meant to represent the decay towards the lin- 
ear torque at large viscosity. We note that Eq. (|151|) im- 
plies, in agreement with the results of the reduced model, 
an overshoot of the edge term of about 20 %, much as the 
one we noticed in the barotropic case (see section IB. 3. ip . 

We note that so far we have assumed that the vorten- 
sity sheet produced at the downstream separatrices has 
the same value as in the unsaturated case. Actually, the 
horseshoe dynamics acts on the entropy in exactly the 
same way as it does with the vortensity in the barotropic 
case, so that the entropy jump at the stagnation point 
can be lower than in the unsaturated case. It is the pur- 
pose of the following section to check the dependence of 
the entropy jump (and therefore of the torque) on the 
diffusivity. Here, in order to assess the correctness of 
Eq. (|150p in our full hydrodynamical runs, we simply 
measure the entropy difference upstream of the U-turns, 
in the saturated state. This difference is the same for all 
runs of set 2, and turns out to be 0.7 times the difference 
in the unperturbed disk. 

Figure Q~7] shows that there is an overall correct agree- 
ment between the results of Eq. (|150p ~ (115ip , and the 
results of numerical simulations. We mention at very 
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low viscosity, the situation is very messy owing to the 
presence of vortices, as described in section 15.3.21 and 
the range of our time average (500 orbits) is too short to 
exhibit a neat trend. We also comment that over the first 
decade of viscosities, the viscous time across the horse- 
shoe region is larger than the duration of the runs, so 
that the corresponding runs have not necessarily reached 
their asymptotic value. 

6.5.4. Dependence on diffusivity 

We now examine how the horseshoe drag edge term 
depends on thermal diffusion. We consider that the 
dependencies of this term on viscosity and on diffusiv- 
ity are separable: its asymptotic value is the product 
of a function of z v exclusively (which we analyzed at 
the previous section) by the load produced downstream 
of the horseshoe U-turns (which depends on the diffu- 
sivity exclusively, as the latter determines the entropy 
jump that subsists upstream of the U-turns). We have 
checked this assumption with an additional set of lower 
resolution runs, not reproduced here, with diffusivity 
K = 3 • 10 - Vf2p. 

We firstly note that the value of the viscosity chosen 
in set 3 {v = 10~ 6 a 2 Q p ) corresponds nearly to the peak 
value in Fig. I17[ so that we make the simplifying assump- 
tion that the horseshoe drag measured in these runs cor- 
responds to the unsaturated horseshoe drag with a lower 
load creation at the downstream separatrices due to a 
lower entropy jump. This yields to the torque expres- 
sion: 



where Ti is the uns aturated horseshoe dra g edge term 
(see e.g. Eq. (96) of IMasset fc Casolj 120091) . and where 
S s , the entropy value at the separatrix upstream of the 
U-turn, is given by Eq. (|109[) . Taking into account the 
cut-off at large diffusivity, Eq. (|152|) eventually yields: 

r « T a min ^1,4^1^4/^ e K « T 1 min(l, \Az]l 2 )e K . 

(153) 

We note that we do not add a blend of the linear torque 
value at large diffusivity, counter to what we have done 
in the analysis of dependence upon viscosity. Although 
we expect to recover the linear torque value at large vis- 
cosity, we expect the disk to behave isothermally in the 
limit of large diffusivity. Since there is neither a vorten- 
sity gradient nor a temperature gradient in the disk, we 
expect the corotation torque to vanish in this limit. 

Figure fTS] shows that there is an overall correct agree- 
ment between Eq. (|153l) and the results of simulations. 
We note that the entropy at very low diffusivity has not 
settled yet to the expected value, as the thermal diffusion 
timescale is longer than the duration of the simulations. 
We also mention that, whereas the correction due to the 
bulk term, given by Eq. (|130j) . was a minute effect in the 
previous section, it is a significant effect here, in partic- 
ular at low thermal diffusion where a significant viscous 
creation of load occurs along the separatrices and gives 
rise to a relatively strong bulk term, as already pointed 
out in sections 15.3.21 and 16.4.11 
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Figure 18. Entropy related torque (edge term) as a function of 
diffusivity. The solid line shows the value given by the simulations, 
corrected from the bulk term of Eq . 11291 . while the dotted line in- 
dicates the prediction of Eq. 11531 . The inset plot shows the ratio 
of the viscosity difference upstream of the U-turns to this difference 
in a n unp erturbed disk (solid line) compared to the expression of 
Eq. H109I) (dashed line). The torque values correspond to time av- 
erages from 800 to 1000 orbits, while the entropy jump is measured 
at r = a ± x s , 4> = ±1 rad, on a time average of the entropy field 
between 700 and 1000 orbits. 

Prior to the full torque expression, which will be given 
in section [8j we briefly present additional results ob- 
tained from a systematic study of the differential Lind- 
blad torque, and we lay down the set of assumptions 
used to provide the final torque expression provided in 
section [5] 

7.1. Diffusivity dependence of the differential Lindblad 

torque 

For the sake of completeness, we have undertaken an 
additional study of the Lindblad torque dependence on 
the diffusivity. At low diffusivity, the acoustic waves 
that constitute the wake essentially behave as adiabatic 
waves, whereas at large diffusivity, they behave isother- 
mally. Since the differential Lindblad torque scales as 
cj 2 , where c 8 is the sound speed, one should recover a 
differential Lindblad torque value 7 times larger in the 
isothermal limit than in the adiabatic limit. We have 
therefore run calculations in which we have a vanish- 
ing vortensity gradient and a vanishing entropy gradient 
(which implies a — 3/2 and / = 1/5) so as to have a 
vanishing corotation torque. These calculations are run 
over a short time interval (typically 10 orbits), which is 
suffici ent for the Lindblad torque to reach a steady state 
value ()Mever-Vernet &: Sicardvl [l987) . We comment that 
we modified the diffusivity procedure for these specific 
calculations, as we had initially linked the amount of 
thermal diffusion to the Laplacian of the entropy, in or- 
der to filter out a possible dependence of the Lindblad 
torque on the diffusivity (see section [2~5]) . In the present 
case we use the standard dependency upon the Laplacian 
of the temperature. In the case of a single acoustic wave 
in a uniform medium at rest, the transition of the phase 
velocity from c s to 7 1 ^ 2 c s depends on the dimensionless 
parameter p = c 2 /(uik), where ui is the wave frequency. 
This parameter represents the ratio of the wave period to 
the diffusion timescale across one wavelength. We there- 
fore expect the transition of the Lindblad torque to occur 
for a critical value of the diffusivity that is 
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(154) 



Saturated torque formula for planetary migration 



23 



obtained by letting p — 1 and adopting u> — Q p for the 
frequency of the waves seen in the matter frame, since 
those are excited at their Lindblad resonances (where the 
wave frequency matches the epicyclic frequency, hence 
the orbital frequency in a Keplerian disk). We note that 
the characteristic length scale in the wake excitation re- 
gion is H, which is also the largest scale of turbulence, so 
that the diffusion approximation should be on the edge of 
its domain of validity, much as for the horseshoe drag, as 
was discussed in section l2~Tl This critical value of the dif- 
fusivity is much larger (by two orders of magnitude) than 
the maximum diffusivity considered in our exploratory 
analysis of the corotation torque. In the Lindblad torque 
runs, we find that indeed the Lindblad torque has a value 
that is roughly the mean of the isothermal and adiabatic 
limits for this value of the diffusivity, and we provide the 
following fit to the dependency that we have obtained: 



the horseshoe width the simple value 



L LR 

-piso 
1 LR 



/(k/k c ,7), 



where f{x,^) is defined by: 

(z/2)V2 + i/ 7 



(x/2)Va + i 
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7.2. Assumptions for the torque expression 
7.2.1. Transition from two dimensions to three dimensions 

The study that we have presented in this paper is 
two dimensional, for reasons of computational cost, so 
that the final torque expression is smoothing depen- 
dent. Useful torque expressions should however be three- 
dimensional, so that we scale our tor que value by the 
three- dimensional estimates given by iMasset & Casolil 
( 2009) . These estimates were obtained by assuming hor- 
izontally layered motion in the coorbital region, and by 
integrating the contribution of all slices, assumed to be 
given, by virtue of Pythagoras'theorem, by two dimen- 
sional estimates in which the smoothing length is the 
layer's altitude. This approach, albeit not genuinely 
three dimensional, is more sophisticated than a two di- 
mensional study with a fixed, reasonable value adopted 
for the smoothing length. 

We note that the saturation of the corotation torque 
depends on the ratio of the libration timescale to the 
diffusion time scale (be it the thermal or viscous diffu- 
sion timescale) . This ratio depends on the altitude in the 
disk, and its estimate requires an estimate of the width of 
the horseshoe region at each altitude. Since most of the 
torque arises from the layers which are located near the 
equatorial plane of the disk, we simply assume that the 
saturation function of the three dimensional torque can 
be correctly described by the saturation of the two di- 
mensional torque in which we adopt the horseshoe width 
at low softening length (i.e. low altitude). While this 
may be marginally true for the vortensity related corota- 
tion torque and more generally for all the bulk terms, it is 
likely a correct procedure for the entropy related torque 
and more generally all the edge terms, the contribution 
of whic h is thought to arise m ainly from the equatorial 
regions (|Masset fc Casolil l2009f) . We therefore adopt for 
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which is close to the maxima l value found by 
iPaardekooper k, Papaloizou ( 2009b). We stress that this 
approximate value is required only to estimate the argu- 
ments of the saturation functions, not the torque mag- 
nitude itself. The uncertainty on these arguments due 
to the uncertainty on the value of x s is much less than 
the uncertainty arising from our ignorance of diffusive 
processes in protoplanetary disks. 

7.3. A note on the cut-off functions 

Accounting for the behaviour of the entropy related 
corotation torque at large diffusion requires to take 
into acco unt the decay of the torque towar ds its lin- 
ear value (jPaardekooper fc Papa loiz ou 20 09a), using the 
cutoff functions evaluated in section [5ATJ We note how- 
ever, as mentioned above, that we do not expect the same 
value in the limit of a large viscosity (where we should 
recover the linear estimate) and in the limit of a large dif- 
fusivity (where we should get rid of the entropy related 
torque, as the fluid becomes isothermal). For the sake 
of simplicity we make the conservative assumption that 
the torque limit vanishes in both cases, motivated by the 
fact that the horseshoe drag edge term should be signif- 
icantly larger than the entropy related lin ear corotation 
torque (jPaardekooper fc Papaloizou 2001). 

Regarding the cutoff of the bulk term, we use a linear 
vortensity related torque value that is 2/3 of the corre- 
sp onding horseshoe drag, by comparing the value given 
bv IMasset fe Casolil (|2009f ) in Eq. (102) to that given by 
Taua kaet alJ (|2002| ). We note that there is no need to 
apply a cut-off at large thermal diffusion to the bulk term 
arising from the entropy- vortensity viscous coupling, as 
this term is negligible in the limit of large thermal diffu- 
sion. 

8. GENERALIZED TORQUE EXPRESSION 

In this section we provide the torque formulae that 
apply to low mass planets embedded in viscous disks 
with thermal diffusion. For the convenience of the reader 
mainly interested in the results rather than in the deriva- 
tion, this section is self-contained. 

All the torque values given below are expressed in 
terms of the canonical torque value r ro f given by: 

r rcf = Z c n 2 p a 4 q 2 h- 2 , (158) 

where the reader is referred to Tab. [1] for the notation. 
The total tidal torque r to t exerted by the disk on the 
planet features two components: the differential Lind- 
blad torque Tlr and the (non-linear) corotation torque 
or horseshoe drag Tqr'- 



r tot = r 



LR 



r cr 



(159) 



The d ifferential Lindblad torque is given bv lTanaka et al.1 
(2003): 



Trcf 



(2.3 + 0.4/3 -0.1a)/ 



(160) 



where k c is defined by Eq. (|154[) and the function / 
is defined by Eq. (|156p . We note that there has been 
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some recent controv ersy on the coefficient of t he tem- 
perature gradient /?. iPaardekooper et all (|2009l ) give in 
their Eq. (14) a coefficient as large as 1.7, which is cor- 
roborated by two dimensional numerical simulations with 
a fixed smoothing length of the potential. We quote 
he re the coefficient 0.4 derived from the data provided 
by iTanaka et al.l ((2002). A definite answer to this issue 
awaits dedicated three dimens ional calculations, such as 
those of lD'Angelo et al.l (|2003|) , with a systematic explo- 
ration of the temperature gradient. Such calculations 
ar e under way, and see m to indicate that the expression 
of Tanaka et al. (2002) is correct, with a high accuracy 
(D'Angelo, priv. comm.). 

The multiplication by the function / in Eq. (jl60[) is 
intended to account for the mild dependence of the Lind- 
blad torque on the diffusivity. When the latter is small 
enough (namely, when the relaxation timescale of tem- 
perature disturbances is larger than the dynamical time) , 
one may substitute the function / by the constant I/7. 

The horseshoe drag features, in general, two terms: a 
bulk term Ths and an edge term r^HS : 



r CR = r 



HS 



Tans- 



(161) 



The bulk term, using Eq. (|129[) a nd applying the norma l- 
ization provided by Eq. (101) of IMasset fe Casolil (|2009D . 
is given by: 



■■7.8z„ 



ref 



VF(z v ) - 2S 



z v F{z v ) - z K F(z K ) 



162) 



-0.62V(1 - £&), 



where z v and z K are defined respectively at Eqs. (j?9"j) 
and (f8Tj) . and where the function F is defined by 
Eq. (|120p . The value of x s , the half width of the horse- 
shoe region, that is needed for the estimate of z v and z K , 
is given by Eq. ([T5?|t . Finally, Eq. ([T62]) also displays the 
bulk cutoff function give n bv Eq. (|70p. and which ca n 
be recast, using Eq. (64) of iBaruteau fc Massetl (|200l . 



e h ~ (l + 30fca^) _1 . 
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The edge horseshoe drag, using Eqs. (|150[) . (|151l) . (I153[) 
and applying the norm alization provided by Eq. (98) of 
IMasset fe Caso li (2009), is given by: 



= —3.35 \Az X J 2 s K 1.8zl^ 2 e„, 



(164) 



where, using respectively Eqs. (|7T|) and (|72j) : 

e, y -[l + (6/iz l/ ) 2 ]- 1 (165) 

e K -(l + 15ftz K )~ 1 (166) 

and where 

X = mm(l,X). (167) 

The set of Eqs. (|159[) to (|167l) provides self-contained 
torque formulae for a low mass planet in a disk with arbi- 
trary viscosity and thermal diffusivity. In the particular 
case in which viscosity and thermal diffusivity remain 
moderate (more precisely when z VjK -C h~ l ), one may 
drop the cutoff coefficients £b,v,n, which are then close to 
unity. 

We comment that the estimate of the bulk term of 
the horseshoe drag, given by Eq. (|162l) . is more accu- 
rate and robust than the estimate of the horseshoe edge 



term, given by Eq. (|164p . which heavily depends on small 
scale properties of the flow in the vicinity of the stag- 
nation point of the hor seshoe region, as was shown by 
IMasset fc Casolil §009) and as was emphasized in sec- 
tion 16.5.11 Also, the factor 3.3 that features in this 
eq uation, which st e ms fr om the conservative estimate 
of IMasset fc Casolil (l2009h . is a minimum value. In any 
case, this factor is large and should be sufficient to halt 
migration at locations where the entropy gradient is neg- 
ative and sufficiently large in absolute value, and where 
the viscosity and thermal diffusivity have adequate val- 
ues to prevent saturation. As a consequence, a minor 
change in the value of this coefficient should have little 
qualitative impact on scenarios of migration. 

9. DISCUSSION 

9.1. On the planetary mass range relevant to this 

analysis 

The torque expressions provided here are strictly valid 
in the regime of small planetary mass, for which the 
torque scales as the square of the planet mass (and the 
horseshoe w idth scales as the sq uare root of the plane- 
tary mass). IMasset et al.l ()2006l ) have shown, by means 
of a two dimensional analysis, that when the planetary 
Hill radius represents some fraction of the disk pres- 
sure length scale (typically such that q/h 3 ~ 0.6 in 
their study), the horseshoe zone width is much larger 
than predicted by the g 1 / 2 scaling, extrapolated from 
lower masses, thereby boos ti ng the corotation torque. 
IPaardekooper fc Papaloizou (2009b) have shown that 
this behavior is due to the wake, which, owing to its prox- 
imity, affects the enthalpy at the horseshoe stagnation 
point, hence it affects the horseshoe width. They found 
this behavior to be smoothing dependent: the mass at 
which the horseshoe drag boost is observed depends on 
the softening length, and is smaller at smaller softening 
length. If one translates this in terms of altitude in a 
three dimensional disk, this suggests that the horseshoe 
region is always over wide in the vicinity of the equator, 
over a vertical length scale that increases when the plane- 
tary mass does. The vertically integrated bulk horseshoe 
drag should then be boosted when the vertical extent of 
the over wide region represents a significant fraction of 
the disk thickness. This expectation yields a relation- 
ship similar to that quoted above, i.e. one expects the 
boost to occur for q/h 3 ~ 0(1), where the dimensionless 
number of the right hand side is still to be determined 
by three dimensional calculations. We comment however 
that the edge term of the horseshoe drag, which is bi- 
ased towar ds low altitude regions o wing to its dependence 
in er 1 (sec Masset & Casoli 2009) should be boosted at 
planetary masses even lower than the bulk term. This 
constitutes another reason why the edge horseshoe drag 
given by Eq. (|164[) should be regarded as a conservative 
estimate, even for planetary masses significantly smaller 
than 0.6ft 3 M». 

9.2. Relevance of a diffusion equation to model the 

entropy evolution 

We have modelled the radiative processes at work in 
the disk by a simple radial diffusion equation on the en- 
tropy. The merits of this simple approach have been 
discussed in section [2~5l The actual processes that con- 
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tribute to the energy balance of fluid elements in the 
disk, howe ver, cannot be r educe d to this simple minded 
approach. Kiev fc Cridal (|2008l) have for instance per- 
formed numerical simulations that take into account, in 
addition to the processes that can be described as a dif- 
fusion, the local processes that affect the energy balance, 
i.e. the viscous heating and the radiative losses through 
the disk photospheres. If one neglects the diffusion and 
takes into account only local processes, then an alternate 
simple model of thermal processes can be achieved by 
imposing a local temperature prescription, together with 
a characteri stic decay time of the temperature distur- 
bances (e.g. lBaruteaul l2008). What matters in assessing 
the magnitude of the entropy-related horseshoe drag is 
the value of the entropy at the separatrices upstream of 
the U-turns. The latter results, in steady state, from a 
balance between libration and relaxation. We note that 
the dimensionless variable z K , which is used to quan- 
tify the corotation torque saturation due to thermal pro- 
cesses, scales with the ratio of the libration time to the 
thermal timescale Tth = x s /k: 



The drift induced by migration over a synodic pe- 
riod tqcl is therefore: 



3 r 

47T T th ' 



(168) 



Although we have not undertaken a systematic explo- 
ration of the torque saturation in the frame work of a 
tempe rature prescription such as the one of iBaruteaul 
(2008), we note that his results are compatible with the 
formulae of Eq. (I164[) in which one adopts the value of 
z K given by Eq. f|168[) . In a more general manner, we 
suggest that Eq. (|168j) be systematically used to assess 
the degree of saturation of the entropy-related corotation 
torque, and that r t h, the relaxation timescale of temper- 
ature disturbances of scale length x s , be determined by 
a prior analysis of the relevant thermal processes. 

9.3. Planet-disk relative drift 

In the calculations performed for the analysis of the 
entropy related horseshoe drag, we have imposed a power 
law for the kinematic viscosity that yields no radial drift 
in the unperturbed disk. We do not expect that a relative 
drift of the disk and planet, given either by the large 
scale accretion of the disk's material onto the primary, 
or migration, or an admixture of both processes, would 
alter our conclusions, for the following reasons: 

• The planet-disk drift turns out to be unimportant 
in the barotropic case. By extension, it should be of 
little importance for the bulk term of the horseshoe 
drag in the general case. 

• The viscous drift of the vortensity sheets at the sep- 
aratrices, over a synodic period, typically am ounts 
to TQi>/a, whereas its viscous spread reads ^/Avtq. 
The ratio of these two quantities is vQ.^ 1 / ax Sl 
which should be small in any realistic disk and for 
any low mass planet for which planetary migration 
is relevant. 

• We can also work out the typical amount of drift 
due to type I migration, at a typical speed d given 
by 

« ~ TTTT- (169) 
n p M p a y ' 
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where //<j = irYiQa 2 /M* is the reduced disk mass. 
Eq. (|170[) shows that the drift is a small fraction 
of x s except in very massive disks, close to their 
gravitational stability limit. Since the drift is a tiny 
fraction of the horseshoe zone width, it may impact 
the degree of saturation only in the cases for which 
the viscous spread of the vortensity sheets is also 
a tiny fraction of the horseshoe width, that is to 
say when the entropy related torque is significantly 
saturated. 

10. SUMMARY 

We have undertaken a study of the corotation torque 
saturation, as a function of viscosity and thermal dif- 
fusion, and we have obtained asymptotic torque expres- 
sions that are given in section[8l A property that is essen- 
tial to our derivation is that the horseshoe drag, whether 
barotropic or not, depends exclusively on the distribu- 
tion of vortensity within the horseshoe region, so that 
our analysis essentially reduces to finding the vortensity 
distribution for a steady flow in the corotating frame. In 
particular, the contribution of the under- or over-dense 
regions that appear within the horseshoe region should 
not be added "manually" to the horseshoe drag estimate. 
Doing so leads to a violation of first principles, as is dis- 
cussed in section 16.5.11 Our analysis shows that the sat- 
urated horseshoe drag can be decomposed into a bulk 
term and an edge term. The latter scales with the en- 
tropy gradient, as in an unsaturated situation, and its 
saturation properties depend both on the viscosity and 
on the thermal diffusivity. The bulk term, however, not 
only depends on the vortensity gradient, as in the unsat- 
urated case, but also on the entropy gradient, as a result 
of the viscous production of vortensity at steep entropy 
gradients. 

We have also contemplated in section 17.11 the transi- 
tion from the isothermal to adiabatic differential Lind- 
blad torque (the latter is a factor 7 than the former), 
and found that it occurs when the thermal timescale is 
of the order of the dynamical timescale. The dependence 
of the entropy related torque on the thermal timescale is 
regulated by a markedly different ratio, namely the ratio 
of the thermal timescale to the libration timescale, which 
gives the degree of saturation. 

At large viscosity or thermal diffusivity, the horseshoe 
drag is found to decay. We have made an attempt to ac- 
count for this decay by the use of ad hoc simple analytical 
fits (cutoff functions). A detailed study of the interme- 
diate regime that settles at large diffusion (be it viscous 
or thermal), and that yields a torque value comprised 
between the horseshoe drag and the linear corotation 
torque could be useful for sub-terrestrial mass objects 
embedded in disks with significant viscosity or short ra- 
diative timescales. One can estimate the value required 
for the viscosity or thermal diffusion to recover the linear 
regim e. W ri ting that e^v k, "C 1, we are left with, using 
Eqs. ([163]) . (fl65|) and (fT66|) : 



a 2 il 



> 0.1g 3/2 /^ 5/2 . 
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For the condition on viscosity, this can be translated into 
a lower limit on the a coefficient of the disk: 

a » 0.1q 3 / 2 h- 9 / 2 . (172) 

Eqs. (|171[) and (|172l) show that the corotation torque 
should approach its linear value for an Earth mass planet 
in a disk with h = 0.04, invaded by magnetorotational 
turbulence (the critical a value of the right hand side of 
Eq. (I172[) being then 10~ 3 ). If the planet has a higher 
mass, if the disk is thinner, or if the planet is in the dead 
zone, the corotation torque should be determined by a 
horseshoe drag analysis. 

Finally, we have underlined that the horseshoe drag 
edge term is potentially extremely strong in a three di- 
mensional case, so that the estimate that we have given 
in Eq. (|164[) should be regarded as conservative. This 
stresses the need for a three dimensional study, that 
would allow a calibration of the unsaturated edge term 
value. 

The numerical simulations performed in this work have 
been run on a 140 core cluster funded by the program 
Origine des Planetes et de la Vie of the French Insti- 
tut National des Sciences de I'Univers. Partial support 
from the COAST project {Computational ASTrophysics) 
of the CEA is also acknowledged. F.M. also wishes to 
thank G. Kocnigsberger for hospitality at the Instituto 
de Ciencias Fisicas of UN AM, Mexico. 
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